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ANALYTICAL AND EXPERIMENTAL PERFORMANCE OF OPTIMAL 
CONTROLLER DESIGNS FOR A SUPERSONIC INLET 
by John R. Zeller, Bruce Lehtinen, Lucille C. Geyser, and Peter G. Batterton 

Lewis Research Center 

SUMMARY 

This report applies the techniques of modern optimal control theory to the design of 
a control system for a supersonic inlet. The inlet control problem was approached as a 
linear stochastic optimal control problem using as the performance index the expected 
frequency of unstarts. The details of the formulation of the stochastic inlet control 
problem are presented. The computational procedures required to obtain optimal con- 
troller designs are discussed, and the analytically predicted performance of controllers 
designed for several different inlet conditions is tabulated. The experimental implemen- 
tation of the optimal controllers is described, and the experimental results obtained in 
the Lewis 10- by 10- Foot Supersonic Wind Tunnel (SWT) are presented. 

The design studies showed that the amplitude- frequency distribution of the distur- 
bance seen by the inlet has a large effect on the performance capabilities of the optimal 
controller. In this study two distinct disturbance spectra were assumed. The results 
show that the more disturbance energy there is at high frequency, the more difficult it 
is to control the inlet. 

The experimental program pointed out certain of the problems involved in imple- 
menting a complex modern optimal controller. Controllers were implemented and eval- 
uated with both analog and digital computer components. Analytically predicted and ex- 
perimental frequency response performance compared quite well. The analog and 
digital computer implementations of a particular optimal controller design showed com- 
parable performance results. Computer routines which were used to implement the 
digital computer version of an optimal controller are included. Recommendations as to 
further activities in using the capabilities of linear stochastic optimal control theory are 
also included. 



INTRODUCTION 


The techniques of modern optimal control theory have been applied to the design of 
a control system for a supersonic inlet. A supersonic inlet is that portion of a super- 
sonic propulsion system which decelerates air from supersonic velocity (relative to the 
aircraft) ahead of the aircraft to subsonic velocity at the entrance to the compressor. 
This deceleration is needed because present compressors require subsonic air to oper- 
ate efficiently. The dynamic head of supersonic air at high Mach numbers may com- 
prise a large percentage of the overall propulsion system compression, and, therefore, 
efficient recovery of the pressure head is a critical part of the supersonic propulsion 
system. For subsonic propulsion systems, however, almost all the compression is 
done by the engine’s compressor. To aid the supersonic inlet in operating at peak effi- 
ciency in the face of varying flight conditions, variable geometry features and associated 
controls are required. 

A typical axisymmetric mixed compression inlet is shown in figure 1 in a normal 
operating configuration. Air at supersonic velocity enters the inlet past a weak oblique 
shock wave. It is compressed supersonically past a minimum area point, or throat, up 
to the terminal normal shock. Thereafter, the flow is subsonic up to the compressor 
face station. 



Figure 1. - Schematic of axisymmetric mixed compression supersonic inlet. 


A stable operating condition for the inlet is one in which the throat Mach number is 
greater than one and the normal shock is downstream of the throat. This is the so-called 
started condition. An upstream or downstream disturbance may, however, cause the 
throat Mach number to drop to one, or it may cause the normal shock to move ahead of 
the throat. When either of these occur, the inlet unstarts and enters an undesirable, 
unstable operating region (called unstart). 

During an unstart a shock wave sweeps out of the throat and a strong shock wave 
forms ahead of the inlet. The result is a large increase in drag and a large decrease in 
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the pressure recovered at the compressor face. In addition, there may exist an oscilla- 
tory flow pattern within the inlet. Such a condition of unstart occurring in flight may not 
only interact with the engine, producing compressor stall and/or combustor flameout, 
but the increased nacelle drag and thrust loss can cause a sudden yawing of the aircraft. 
Control is required to maintain throat Mach number and terminal shock position within 
acceptable limits while maintaining efficient inlet operation. 

Basic control devices are bypass doors and a variable centerbody. Opening the 
bypass doors allows air to be dumped overboard, causing the shock to move downstream 
away from the throat region. The movable centerbody varies the throat area, thereby 
varying the throat Mach number. A proper combination of these two control variables 
is used to ensure stable (started) inlet operation in the face of upstream and downstream 
disturbances. 

Inlet control systems (refs. 1 and 2) have been designed to minimize system re- 
sponse to deterministic disturbances. Designs were obtained using frequency domain 
techniques. In reference 3, Barry conducted a design study based on an explicit de- 
scription of inlet disturbances. The disturbance treated was atmospheric turbulence 
described by experimentally determined power spectral densities and probability dis- 
tributions. The criterion used for evaluating inlet controls was the expected frequency 
of inlet unstart. 

The control system to be discussed in this report has been designed to minimize 
the unstarts that would be initiated by a downstream (engine compressor face) airflow 
disturbance. This approach is an extension of the work of Barry. Initial work in this 
area has been presented in references 4 and 5. The inlet control problem was ap- 
proached as a linear stochastic optimal control problem using, as the performance 
index, the expected frequency of unstarts. References 4 and 5 document the theoretical 
basis and computational procedures required in designing and analytically evaluating 
modern optimal inlet controllers. 

The techniques of modern optimal control theory as applied to inlet control design 
are being investigated for several reasons. First, the modern approach provides a 
rigorous solution technique for optimizing a control design to some specific performance 
criterion. Second, the resulting control design will be stable. Stability is not neces- 
sarily assured when using conventional techniques. Third, the approach is general 
enough that system constraints can be included in the performance criterion. For ex- 
ample, in the inlet problem, limitations on bypass door position, velocity, and acceler- 
ation can be taken into account by a proper formulation of the criterion. Fourth, noisy 
measurements as well as random disturbances fit quite well into the modern optimal 
control formulation. Finally, the theory is such that it can handle the multiple- input - 
multiple- output control problem. The inlet control design, although it is not so con- 
sidered in this report, can be expanded to a multiple input- output problem. Such an 
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approach could sense additional pressures and pressure ratios throughout the inlet duct 
and control with centerbody as well as bypass doors. 

The analytical inlet model used for the controls analysis of reference 4 was a sim- 
plified representation of an actual experimental variable geometry mixed compression 
supersonic inlet under evaluation at Lewis Research Center (refs. 6 to 8). Controller 
performance was evaluated analytically for wide spectrum (white) stochastic disturb- 
ances at two different levels of measurement noise on the sensed output variable. This 
report, however, expands the inlet model to include additional aspects such as (1) the 
response limitations of the actuators for the control input (overboard bypass doors) and 
(2) downstream airflow disturbances which have nonwhite (colored) power spectral den- 
sities. These considerations are required when the controller designs are to be im- 
plemented and evaluated experimentally on an inlet operating in a supersonic environ- 
ment. 

The rigorous solution techniques of modern optimal control design generate a con- 
troller configuration which in most cases is considerably more complex than one ob- 
tained by classical cut-and-try frequency domain techniques. It is the purpose of this 
report to discuss the procedures involved in determining the optimal design and then im- 
plementing with hardware the complex modern controller configuration. In addition, in 
determining the linear state- space inlet representation, several approximations to the 
real nonlinear distributed parameter inlet model have to be made. This report, there- 
fore, uses results from the experimental operation of the controllers to determine the 
adequacy of these approximations. 

The information is presented in two parts. First, the details of the formulation of 
the stochastic inlet control problem are discussed and documented. Along with this is a 
description of the computational procedures required to arrive at the optimal controller 
designs. Finally, a tabulation of the analytical results of controllers for several differ- 
ent inlet conditions is presented. In the second part, the details of the hardware imple- 
mentations of controllers are described. This is followed by a presentation of the ex- 
perimental results obtained in the Lewis 10- by 10- Foot Supersonic Wind Tunnel (SWT), 
as well as a comparison of these results with the frequency responses as predicted ana- 
lytically. Finally, some recommendations are presented as to further efforts that would 
enhance this initial endeavor at applying optimal controller design to supersonic inlets. 


CONTROLLER DESIGN AND ANALYTICAL PERFORMANCE 
General Solution Technique 

As stated earlier, the design techniques of modern optimal control theory have been 
applied to the design of a control system for a supersonic inlet. The purpose of the inlet 
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control system considered here is to minimize the expected frequency of inlet unstarts to 
a random downstream airflow disturbance. Figure 2 is a block diagram of the general 
configuration chosen for this study. 

As shown in figure 2, there are three distinct transfer functions, two defining the 
inlet and one the downstream (compressor face) disturbance. They are (1) the dynamics 
of the subsonic duct designated as (2) the dynamics of the bypass doors to 

be used for control designated as Ggp D (s), and (3) the transfer function G^g(s) (noise 
shaping network) which models the dynamics of the airflow disturbance to the duct. 

Each of these is discussed in detail in later sections of this report. 
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Figure 2. - Block diagram of typical inlet configuration. 


The Gaussian compressor face disturbance w d , shown in figure 2, is modeled as a 
white Gaussian airflow disturbance w being operated on by the transfer function G Ng (s). 
The control input u operates the bypass doors and results in a corrective control air- 
flow w . A measurement z of terminal normal shock position y_ is measured 

C © 

through a noisy channel with measurement noise v . The measurement noise is 

y s 

assumed to be white Gaussian. 

For the inlet control design the following performance index was chosen to be mini- 
mized: 


J = A + k ctJ 


where 


( 1 ) 


A 


2n 



( 2 ) 
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and 

X expected frequency of inlet unstarts 

k positive weighting factor 

or mean- square value of control input 
o 

o? mean- square value of shock velocity 

y s 

o 

a mean square value of shock position 
y s 

a shock position tolerance (distance between undisturbed shock position and inlet 
throat, see fig. 1) 

The cost J was selected so that the control must minimize unstarts X while limiting 

2 

the amount of bypass door control effort a u needed to do so. (All symbols are defined 
in appendix A. ) 

The X relation (eq. (2)) gives the expected frequency with which the Gaussian ran- 
dom variable y„ exceeds the level a in the positive direction. The derivation of 
equation (2) can be found for instance in reference 9. The weighting factor k for a* 
is selected to penalize the control variable so that the level of control effort will not 
exceed that which is available. (Selection of the control effort weighting factor k is 
discussed in a later section. ) In order to use X of equation (2) for this control design, 
the following assumptions must be made: (1) the inlet disturbances are Gaussian, (2) the 
inlet dynamics are linear, and (3) the controller is restricted to being linear and time 
invariant. 

The approach taken in the designs being presented in this report uses the techniques 
of linear stochastic optimal control and estimation theory. This solution involves min- 
imizing a quadratic type of performance index. It should be noted that the performance 
index of equation (1) is not quadratic because of X. A linear optimal control law can, 
however, be determined by employing a technique termed the quadratic equivalence prin- 
ciple (ref. 10). This technique is used in this report. Since it has been previously de- 
scribed and used in reference 5, it is not repeated here. A summary of the type of 
control system which results is presented in the following section. 


Linear Stochastic Optimal Control and Estimation Solution 

A linear time invariant system can be described in state variable form as 

x = Ax + Bu + Dw (3) 
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where x is an n x 1 state vector, u is a cxl control vector, and w is a d x 1 plant 
disturbance vector. An l X 1 output vector y is defined as 

y = Cx (4) 

and an m x 1 measurement vector is defined as 

z = Hx + v (5) 

where v is an m x 1 measurement noise vector. Both w and v are white zero mean 
Gaussian and uncorrelated with each other. Quantities A, B, C, D, and H are matri- 
ces of appropriate dimensions. 

In solving the control and estimation problem (using the approach of ref. 11) for a 
quadratic performance index, the following equations result. The feedback control law 
is defined as 


u = -K c x (6) 

where x is the optimal estimate of the state vector x and is generated with a Kalman 
filter described by 


x = Ax + Bu + K 0 (z - Hx) (7) 

The computation details for the constant matrices K c and K g are given in reference 4. 

The block diagram in figure 3 illustrates the solution to the optimal control and es- 
timation problem showing the state estimator and state estimator feedback. The state 
estimator (Kalman filter, eq. (7)) is basically a model of the plant driven by control u 
and measurement z. Signal z is compared with the estimated measurement z to 
form a term which is the error in the estimate of the measurement. This error is then 
multiplied by Kalman filter gains K 0 and added back into the filter as a ’’correction" 
term . The filter output x is weighted by the control gains K c to form the optimal 
control vector u. The portion of the system with measurement z as the input and con- 
trol u as the output is defined as the optimal controller. 

The following sections discuss the details of the design procedure as applied to the 
inlet problem. 
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Figure 3. - Combined optimal regulator - state estimator. 

Linear Continuous Time-Invariant Model Formulation 

Inlet transfer functions . - The experimental mixed compression inlet, for which an 
optimal controller 1ms been designed, has been the subject of evaluation in several past 
programs at Lewis Research Center, During these programs, the dynamic relations 
between a downstream (compressor face) disturbance and specific measurable variables 
throughout the inlet have been determined. These relations have been obtained through 
a frequency response testing method (refs. 6 and 8). In appendix B is a brief descrip- 
tion of the method and how it was used in evaluating the inlet open loop frequency re- 
sponse performance. Appendix C contains the frequency response data obtained to de- 
scribe the performance of the inlet. Also in appendix C is a complete tabulation of the 

transfer functions which have been curve fit to the experimental data. These transfer 

/ - TjS\ 

functions involve transportation lags or pure dead-time terms (e a J. This is to be 
expected considering the distributed nature of the inlet duct. For comparison purposes, 
appendix C contains the frequency responses of the transfer function approximations to 
the experimental data. 

For the inlet controls program being documented in this report, two specific meas- 
urements of shock location have been considered. One configuration uses a sensor 
which provides a stepwise continuous indication of actual inlet shock position. This 
type of sensing of the normal shock position has been accomplished in previous research 
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programs (refs. 12 and 13). A brief description of the technique is also included in 
appendix B. The control using this measurement of actual shock position is designated 
as the shock position feedback or SPF system. It is described by the block diagram in 
figure 4. The second shock measurement configuration uses a static pressure down- 
stream of the throat to indicate the position of the inlet normal shock. This is a more 
conventional way of obtaining an indication of shock position. For this second configura- 
tion, it was assumed that only the throat exit static pressure p te was measurable for 
purposes of control. Thus, the inlet is uncontrolled or open loop to the actual shock 
position location. This configuration, which is shown in the block diagram in figure 5, 
shall be designated as the throat exit feedback or TEF system. 

Noise shaping 



Figure 4. - Block diagram of shock position feedback (SPF) system of 
inlet control. 


Noise shaping 


i 



Figure 5. - Block diagram of throat exit feedback (TEF) system of inlet control. 
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For the SPF system (fig. 4), the transfer function relation of shock position y_ in 
response to a downstream (compressor face) airflow disturbance is given as 


G INLEt( S ) 


ifi or / s i\ -4.0x10 s 

16 - 25 Ur 

w * (s) /i + VjL + m +1 \ 

\80 /( 365 2 365 ) 


cm 

kg/sec 


(Cl) 


Since our intent is to develop a finite- order state variable formulation of the inlet, the 
dead- time term of equation (Cl), which has an infinite number of poles, must be modi- 
fied. A finite- order approximation for the dead time was obtained using a Pade approx- 
imation. The transfer function for the SPF system inlet model can be written as 


G INLET^ 


Yg(s) 

Wj(s) 


16. 25 + 1)(1 - 2xl0' 3 s + 1. 6x10" 6 s 2 - 0. 533xl0" 9 s 3 ) 

\210 l 

(— + lV-4 + + 1|(1 + 2x10" 3 s + 1. 6x10" 6 s 2 + 0. 533xl0" 9 s 3 ) 

\80 /\365 2 365 / 


cm (8) 
kg/sec 


using a third-order Pade. 

For the TEF system (fig. 5), the two transfer functions are 



G SHOCK <s) 


r 8 < s > 

P te (s) 


10. 68e" 2-5x10 3 ® 



cm 

N/cm 2 


(C3) 
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Again, as in the case of the SPF system, Pade approximations to the delay terms 
are used. Since the delay terms in equations (C2) and (C3) are both of shorter duration 
than the total duct delay of equation (Cl), it was determined that first-order Pade ap- 
proximations provided sufficient accuracy. The resulting transfer functions are 


G DUCT^ ~ 


Pte (5) 

w i (s) 


1. 52 /jL + lV-®_ + l)(l - 7. 5X10' 4 S) 
\210 /\500 / 

/_§_ + + l\(l + 7. 5x10" 4 s) 

\80 /\365 2 365 / 


N/cm 2 

kg/sec 


(9) 


^ _ y s (s) _ 10.68(1 - 1.25xl0“ 3 s) 

g shock' s ' ~~T17 ) -7~ s \ 

p te' s ' /-§- + 11(1 + 1. 25x10 d s) 

\500 / 

Bypass door transfer function . - The mechanism used as the control input for the 
mixed- compression inlet under investigation are overboard bypass doors. These are 
fast-acting high-performance devices and are discussed in reference 14. Frequency 
response data from reference 14 are displayed in figure 6. As can be seen, the dynam- 
ics are not linear in that the performance varies as a function of the disturbance ampli- 
tude. Previous tests, however, have shown that a disturbance equivalent to the 14- 
percent level of door movement moves the shock position over a range quite adequate for 


cm 


N/cm' 


( 10 ) 



Figure 6. - Frequency response of Inlet control bypass doors for three disturbance 
amplitudes. 
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investigation of inlet controller concepts. At this level of bypass airflow, the transfer 
function of equation (11) adequately describes the bypass door performance: 


g bpd( s ) “ 


w e (s) 

u(s) 


0.9435 

s 2 ! 2(0. 5)s [ x 
628 2 628 


kg/ sec 
V 


(ID 


Disturbance noise assumptions . - It has been pointed out earlier in this report that 
the inlet controllers were designed to minimize the expected frequency of unstarts to 
downstream disturbances. This is a statistical performance criteria and involves the 
mean- square value of shock position and shock velocity (ref. 4). Thus, a statistical 
description of the downstream disturbances is required. The linear stochastic optimal 
control theory formulation demands that the disturbance w be white Gaussian noise with 
zero mean. For the inlet problem, the disturbance w^ was not white. To model the 
spectrum of w^, transfer functions were selected to shape a white noise input w. The 
presence of the required shaping transfer functions is shown by the dotted blocks in fig- 
ures 4 and 5. 

At the time of this program, no data were available to define the specific shape of 
the w^ spectrum. Therefore, two different disturbance spectra were assumed. Their 
asymptotic representations are shown in figure 7 . The spectra were selected to allow 



Figure 7. - Power spectral density of disturbance w d (asymptotic representa- 
tion). 
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the comparison of resulting optimal controller designs over as wide a range as seemed 
reasonable. To model the spectra shown in figure 7 as well as to allow for some flexi- 
bility in modeling other spectra in the future, the following generalized transfer function 
was used: 


g ns( s ) -■ 


w d (s) 

w(s) 


+ 1 


, &< 


5 2 (« 1 + a 3 ) 


s + 1 


^ 1^3 


o^g 


(12) 


For the spectra of cases A and B in figure 7, the a parameter values are shown in 
table I. To serve as a basis of comparison, it was decided that the mean- square value 
of the disturbance airflow would be the same regardless of the frequency spectrum 
selected. This was accomplished by modifying the power spectral density level of the 
white noise input w in accordance with the particular frequency spectrum selected. For 
the case A and case B spectra, the white noise input power spectral densities PSD(w) 
required to provide a constant mean-square airflow equal to 0.0282 (kilogram 

2 

per second) are included in table I. 


TABLE I. - FREQUENCY SPECTRUM PARAMETERS 


FOR DISTURBANCE AIRFLOWS 


Parameters 

Case 


A 

B 

Noise shaping transfer function 
parameters, rad/sec: 



“i 

0.1 

0.1 

“2 

5000 

10 

"3 

5000 

2000 

White noise input power spectral 
density, PSD(w), (kg/sec)^ 
(rad/sec) 

0.554 

0. 188 

Mean- square value of disturbance 

airflow, , (kg/sec) 2 

d 

0. 0282 

0. 0282 
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Measurement noise descriptions. - The signals of actual shock position y„ and 
throat exit static pressure p te are used as the output measurements for the SPF and 
TEF systems, respectively. It has been assumed that these measurements are cor- 
rupted by specific levels of additive white Gaussian zero mean noise. This assumption 
is made at this time, since no spectral information of these measurement signals is 
available. In addition, as is discussed briefly in appendix B and in detail in refer- 
ence 12, the shock position sensor generates a stepwise continuous representation of the 
location of the normal shock. No attempt has been made to include the quantization 
error of the sensor in the analysis. The noise levels assumed for this design study are 

PSD ( V y ) = 3 ‘ 22xl0 ~ 4 cm 2 y/ rad/sec 
PSD^v ^ = 2.38 (n/cm j rad/sec 

State Space Model Formulation 

The transfer functions representing the inlet dynamics for the two configurations 
can now be formulated. The inlet frequency domain representations were transformed 
to the state variable form (eqs. (3) to (5)) by using the phase variable transformation 
(ref. 15). Tables n and III are the resulting numerical values for the matrices and 
and vectors for the SPF and TEF systems, respectively. The blocked-in sections come 
directly from the transfer functions indicated on the right. The nonblocked-in elements 
are the coupling between the transfer functions. The values of a^, and for the 
case A and case B disturbance noise spectra are presented in table I. 

It should be noted that both the SPF and TEF systems when put into the state space 
formulation are described by ten first-order differential equations (eq. (3)). Thus, the 
optimal inlet controller (eqs. (6) and (7)) consists of ten gains (K g ) defining a Kalman 
filter which generates ten state estimates (x), which are weighted by ten feedback values 
(K c ). 

With the state-space models of tables n and in the optimal controller design ap- 
proach described earlier (documented in refs. 4 and 5) can now be undertaken. The 
computational details of this procedure are described in the next section. 
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TABLE III. - THROAT EXIT FEEDBACK (TEF) SYSTEM STATE -SPACE MODEL MATRICES 



Blanks are all zeros 

















Computational Design Procedures 


It was desired that the optimal inlet controller be in the form of a combined 
controller-estimator as shown in figure 3. The design procedure is identical to that 
described in detail in reference 4. Therefore, the steps required to determine the esti- 
mator gains (K g , eq. (7)) and optimal controller feedback gains (K c , eq. (6)) are only 
summarized in this report. The estimator gains and the covariance matrix of the esti- 
mation error can be determined by solving an appropriate steady-state matrix Riccati 
equation using the inlet model and noise PSD's. To obtain the control gains K , a state 

V 

regulator problem must be solved. The regulator has the task of minimizing the devia- 
tion of the appropriate states so as to accomplish the minimization of 

J = A + kaj (1) 


when subjected to well defined compressor face airflow disturbances. As stated earlier, 
J involves the expected frequency of unstarts A as well as the effort required to re- 
duce the expected frequency of these unstarts. As was pointed out, A is not a quad- 
ratic term. Thus, to use the results of linear stochastic optimal control theory and 
obtain a linear feedback solution, the quadratic equivalence principle is used. This 
principle is briefly outlined here. 

Consider a general quadratic index in the variables of equation (1): 


J 4i’ ^ s +w 4 s +w 2\) 


The differential of <T is simply 

eq 


dJ eq = + W l da ? s + W 2 d 4 


(13) 


(14a) 


Similarly, the differential of J in equation (1) can be written as 


dJ 



3J /dcr 


dJ/dcrf 


da“ 


dJ/dvl 

3J /do? 

y s 



(14b) 
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it can be seen by comparing equations (14a) and (14b) that dJ = 0, which indicates J 
has been minimized using the same gains K that minimized J . Explicit conditions 
that must be satisfied for J to be minimized can be obtained from the aforementioned 
expressions for W ^ and W 2 by substituting for the required partial derivatives : 



The computational technique for finding the minimum J is outlined in figure 8. There 
will be a minimum J for each value of the control weighting k. The procedure shown 
in figure 8 is repeated for each value of k. 

To determine the minimum cost (J), trial pairs of and Wg (designated as W| 
and W| ) are used as inputs to the optimal regulator portion of the solution. The feed- 
back gains K c are determined by solving a steady- state matrix Riccati equation. Using 
these gains and the covariance matrix of the estimation error, the steady- state state- 
covariance matrix equation is solved to determine mean- square values of the states. 
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Figure 8. - Computation procedure flow chart. 


(This involves solving a Lyapunov equation. ) The mean-square state information is 

used to determine of , of , and of values. These are used to compute J and X as 

y s y s u 

well as and W 2 . When the values of and W 2 computed by equations (15) are 
equal to the trial W| and W| values, then equivalence is achieved and the cost J is 
at a minimum value for that value of control weighting k. 

A search routine on and W| could have been used to find the minimum cost. 
However, it was decided that selected pairs of W| and W| be used, which would en- 
compass the field of possible values, and that both the optimum and nonoptimum solu- 
tions would be printed. This list was then searched manually to find the optimum solu- 
tions. The search was done to see the full deviations of the nonoptimum solutions from 
the optimum. 
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Analytical Design Results and Discussion 


A family of optimal controllers has been designed for each of the four systems dis- 
cussed earlier. These controllers are as follows: 

(1) SPF system, case A disturbance 

(2) TEF system, case A disturbance 

(3) SPF system, case B disturbance 

(4) TEF system, case B disturbance 

Analytical results of the design procedure are presented and discussed in this section. 

If the inlet were left open loop and the undisturbed steady- state position of the shock 
(o;) set by a fixed bypass door opening, then the unstart performance shown in figure 9 
would result. Note that the ordinate is the inverse of the frequency of unstarts. The 
mean time between unstarts X~ * in hours should be a more understandable numerical 
quantity for the reader . 
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Figure 9 is a log- linear plot of X against a for both case A and case B distur- 


bances. It should be remembered that 


w ■ 


is the same for both cases. The straight 


1 2 

lines are due to log (X ) being a linear function of a . This can be seen by examining 
the equation which defines X (eq. (2)). The a = 0 intercepts for the lines are given by 
(see eq. (2)) 


2r Vi/i 

This is the inverse of what is termed the '’zero crossing frequency. " This intercept is 

2 2 

smaller for case B than for case A because of the relative magnitudes of a and ct- . 

/ 2 V 1 Ys ys 

The slope of a line is (2cr ) . This quantity is smaller for case A; thus, mean time 

\ y s/ / 2 Y 1 

between unstarts for case A is less sensitive to a than in case B. The reason (2a y j 

is smaller for case A is that in case A most of the energy is concentrated at low frequen- 
cies where the disturbance energy is not greatly attenuated by the inlet duct. Converse- 
ly, for case B, more disturbance energy is present at high frequencies where the inlet 
attenuation is large; hence, the resulting mean- square shock position is less than for 
case A. 

Figures 10 to 13 represent the inlet unstart performance for each of the inlet prob- 
lems being evaluated, j On each of these figures the ordinate is the time between unstarts 
X" 1 and covers the sable range as thai; of the shaded area of the open- loop performance 
shown in figure 9. 

The abscissa for these four figures is a , the rms control effort in volts. This 
factor was part of the performance index of equation (1). Each value of a u corresponds 
to a different value of control weighting k and thus a different set of feedback gains K c . 
In looking at the curves for any fixed value of shock setting a, the time between unstarts 

increases as the amount of control effort a increases. Also, for a fixed control effort 

-1 U 
but an increased a, X is greater. 

In figure 10 for SPF case A there is a sharp increase in X" 1 in the area of a u = 

0. 175 volt. Beyond 0. 20 volt, no significant gain in performance can be accomplished. 
This type of control performance is also seen in figure 11 for TEF case A which is the 
same disturbance case. When comparing figures 10 and 11 it can be seen that SPF sys- 
tem can use lower a settings to accomplish the same unstart performance. Thus, for 
the case A disturbance, control can better be accomplished by directly sensing the out- 
put y g even though an additional measurement lag is incurred in doing so. This is due 
to the relatively high measurement noise level present on the p^ e measurement signal. 

Also shown in figure 10 is a notation which indicates the shock setting that would be 
required for an open- loop system to yield 100 hours between unstart. At this setting of 
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Shock position 










15 centimeters, inlet overall performance (pressure recovery and distortion) would be 
considerably worse than with the setting of 0.89 centimeter possible with closed- loop 
control. 

Figure 12 and 13 are unstart performances for the two feedback configurations 
(SPF and TEF) for the case B disturbance (high frequency content). The three a 
settings are the same for each configuration. When comparing these two figures it can 
be seen that there is a benefit from sensing the throat exit pressure p^ g instead of 
shock position y . This signal is closer to the disturbance than shock position y , 
and, for the high-frequency (case B) disturbances, it shows unstart improvement. Even 
though the p te signal is highly corrupted by measurement noise, the data of figure 13 
show that the closeness of p te to the higher frequency disturbance yields performance 
benefits. The need for less phase lag between disturbance and measurement seems to 
outweigh the measurement uncertainty where higher frequency disturbances are con- 
cerned. 


Experimental Controller Selection 

In the preceding sections the techniques for finding an optimal inlet controller for 
a nonquadratic performance index were presented. These techniques were applied to 
the design of optimal controllers for the 40/60 inlet, which were then evaluated in the 
SWT. Selection of the feedback gains (K ) and estimator gains (K ) for three of the ,! 

v C ) 

four plant/noise configurations was made. The three configurations selected are 
SPF system, case A disturbance (fig. 10) 

TEF system, case A disturbance (fig. 11) 

TEF system, case B disturbance (fig. 13) 

For each of the three configurations one specific set of optimal controller gains 
corresponding to a specific value of rms control effort ct u was selected. The values 
of a u at which the controller gains were selected are indicated in figures 10, 11, and 
13. The actual selection was carried out in the following manner. 

The physical variable used in selecting the control gains K c was the control bypass 
airflow w . This variable has a well defined maximum value, determined by the maxi- 

V 

m um opening of the bypass doors. The airflow w c is related to the bypass door actu- 
ator input u by the transfer function of equation (11). For each optimum controller, 

the rms value of control airflow or was computed. The controller gains selected for 

w c 

experimental evaluation were those for which the resultant value of a was equal to 

w c 

0. 168 kilogram per second. This value is equal to 10 percent of the rms bypass door 
flow capacity when operating about midposition. 

It should be pointed out that the transfer function of equation (11) indicates the con- 
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trol airflow w c is related to the input u as a function of frequency. Also, the inlet 
disturbance energy has some frequency distribution which causes some type of frequency 
distribution on the control signal u. 

Therefore, even though the three experimental optimal controllers were selected 

for the same value of rms control airflow cr , the rms control efforts a required to 

w c u 

produce this fixed value of rms control airflow are different. This is indicated by the 
different values of a u appearing at the selection points in figures 10, 11, and 13. The 
vertical lines on these figures can be used to determine the unstart performance for the 
selected control designs. 

The experimental results obtained using these selected controllers are presented 
and discussed in the EXPERIMENTAL CONTROLLER PERFORMANCE section. 

EXPERIMENTAL CONTROLLER PERFORMANCE 
Analog (Continuous) Controller 

The state estimator - optimal controller configuration discussed in the CONTROL- 
LER DESIGN AND ANALYTICAL PERFORMANCE section and shown schematically in 
figure 3 is described by the vector-matrix equations 

k = Ax + K e (z - Hx) + Bu (16) 

and 


u = -K c x (17) 

These equations can be implemented directly by using analog computer components. Ap- 
pendix B gives a brief description of the actual computer equipment employed in the ex- 
perimental facility. As discussed earlier, three different optimal controller designs 
were experimentally evaluated. These designs involved two different measured vari- 
ables described earlier as the SPF and TEF systems. Figures 4 and 5 show the general 
block diagrams for these two different configurations. 

As can be seen from equations (16) and (17) and figure 3, the optimal controllers 
lend themselves quite naturally to hardware implementation with an electronic analog 
computer. The only difficulty involved in finalizing the analog hardware involved scaling 
the large values of estimator gains and control gains resulting from the design procedure. 
Most of the scaling problems were eliminated by using very fast integrating rates on 
each of the integrators. In the analog circuit, the outputs of various integrators were 
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the system estimated states x(t). In this particular problem, the presence of transfer 
function zeros and the use of a Pade approximation for the dead time caused the state 
variables to differ from any actual system variables. Thus, the level or magnitude of 
state estimates during experimental system operation could not be easily predicted. 

This made analog scaling to prevent amplifier overloads very difficult. To resolve this 
problem, the experimental analog controllers were operated first with a linear analog 
simulation of the inlet system transfer functions. Scaling the estimated states was then 
adjusted to allow all the amplifiers to operate at satisfactory levels under the worse case 
levels of disturbance inputs. The results of operating the three designs with the experi- 
mental inlet in the wind tunnel (SWT) are presented in a later section. 


Digital Computer (Discrete) Controller 

Since a digital computer was already available in the SWT facility for a companion 
experimental controls program (ref. 16), it was decided to implement optimal control 
laws with a discrete controller. A brief description of the digital equipment is included 
in appendix B and a detailed description is in reference 17. 

Presently the optimal inlet control is formulated and designed as a continuous con- 
troller. Equations (16) and (17) describing the optimal controller are linear, time- 
invariant, differential, and algebraic equations. Two possible approaches for designing 
optimal inlet controllers are available. One method involves transforming the inlet 
open- loop differential equations into discrete-time (difference) equations. Then a com- 
plete optimal control system can be designed in discrete time . Such an approach was 
not used, since for the inlet it would have required a new formulation of the optimal con- 
trol solution as well as the development of new computer routines. 

The other method for obtaining a digital computer control law involves approximat- 
ing the continuous control law of equations (16) and (17) by difference equations. The 
performance of a system using a digital computer to implement these equations can be 
made equivalent to that possible with the continuous controller. This second method is 
the one selected for the program discussed in this report. The block diagram in fig- 
ure 14 shows the manner in which the complete digital computer control system was im- 
plemented. 

First, the appropriate measured output is sampled and converted to a digital equi- 
valent upon which the computer can operate. During the uniform sampling interval T 
the computer algorithm is exercised and an optimal controller output u is obtained. 
This then is input to the control doors at the next sample time and it is held fixed for 
the duration of the sample period T. A sampling period of 1 millisecond was conser- 
vatively selected using the closed- loop stability criteria discussed in appendix D. 
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Figure 14. - Digital computer inlet control system block diagram. 


Systems in which input and output controller information is sampled are called 
sampled- data control systems. Various techniques for analyzing such systems can be 
found in the literature. Both frequency domain approaches using the z- transform 
(refs. 18 and 19) and time domain approaches (refs. 15 and 20) have been used. Since 
the continuous formulation for the inlet problem is in the time domain, a time domain 
discrete formulation was obtained by using the state transition matrix (refs. 15 and 20). 

For this particular control problem with ten estimated states which are not closely 
related to distinct physical variables, certain numerical programming problems were 
encountered. Techniques used to overcome these problems and arrive at an acceptable 
control algorithm are not included in this section since such details are not essential to 
a discussion of the experimental results . However, appendix D is included to discuss 
in detail the techniques involved in arriving at a practical computer control law. 

The experimental performance of the discrete controller is presented in the next 
section along with the analog or continuous controller results. 


Experimental Results and Discussion 

In the experimental program, the system was subjected to sinusoidal airflow dis- 
turbances at the downstream end of the inlet as described in appendix B. Thus, all the 
performance data to be presented are in the form of closed- loop frequency responses of 


\ 


A 
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shock position to an airflow disturbance. The results show how well the controlled sys- 
tem regulates against sinusoidal disturbances of fixed amplitude at different distinct 
frequencies. 

The test program would best have been run with random airflow disturbances as 
described by the spectral densities shown in figure 7. This was not done since the dis- 
turbance devices (bypass doors) would not have been capable of accurately duplicating 
these spectral densities. Also, to measure the frequency of unstarts to a random dis- 
turbance would have required considerable running time. This is impractical in the 
SWT. 

Performing frequency response tests with fixed amplitude sinusoidal signal inputs is 
the technique used in past programs to evaluate controller experimental closed- loop 
performance. This is a direct way to look at linear time- invariant systems. 

The experimental data presented in figures 15 to 23 consist of closed- loop frequency 
responses for the three controller designs discussed in the CONTROLLER DESIGN AND 
ANALYTICAL PERFORMANCE section. Also included are experimental open- loop fre- 
quency responses to evaluate the different controllers. For each controller, compari- 
sons are made between the experimental and analytically predicted closed- loop frequen- 
cy response performances. Analytical predictions of closed- loop performance are 
obtained as follows. Open- loop models of the inlet are defined by the finite- order trans- 
fer functions determined in appendix C and represented in the time domain by the matri- 
ces of tables II and III. A closed- loop system transfer function is derived and the fre- 
quency response is evaluated by using the open- loop models and the appropriate optimal 
gains K c and state estimator. Also, where possible, comparison is made between the 
analog and digital computer implementations of the control laws. 

Only frequency response magnitudes, not phase angles, are presented. All magni- 
tude data are normalized to the open- loop magnitude at 1 hertz. 

SPF system frequency response . - Figure 15 is a frequency response plot of the 
SPF case A controller design implemented with analog computer components. The un- 
controlled or open- loop response of inlet shock position to a downstream airflow disturb- 
ance is also shown. It can be seen that control attenuates shock motion by a factor of 
at least 10:1 at frequencies of 0.5 hertz or less. However, as the disturbance frequency 
increases, the controller fails to attenuate the disturbance as well. In fact, from about 
6 to 20 hertz it even amplifies the disturbance somewhat. Beyond this frequency, the 
shock position behaves as if the system were open loop. The case A disturbance, for 
which this particular control was designed, contains the majority of its disturbance 
energy at low frequencies. Thus, a large closed- loop attenuation is produced at the 
lower frequencies. The control design assumes what little disturbance energy exists 
at high frequency is sufficiently attenuated by the inlet duct dynamics; thus, the closed 
loop follows the open loop in this area. The slight amplification of shock motion from 
about 6 to 20 hertz probably does not significantly increase unstart frequency, since 
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disturbance energy in this band is small. 

Figure 16 is a comparison of the analog SPF case A control experimental perform- 
ance with the analytically predicted closed- loop performance. It can be seen that the 
analytical predictions show more attenuation at low frequency than the experimental 
values. It was found during the SWT tests that the shock position gain was about two- 
thirds the value used in modeling the plant (numerical values of table n) and designing 
the controller. This gain discrepancy is the most probable cause for the difference be- 
tween the analytical and experimental performances especially at low frequency. 

Figure 17 is a comparison between the experimental performance of the analog (con- 
tinuous) and a digital computer (discrete) implementation of the SPF case A control de- 
sign. The two implementations are, in general, quite similar except for some lack of 
low frequency disturbance attenuation with the digital version. 

TEF system frequency response . - Figure 18 shows the experimental closed- loop 
frequency response of shock position to a disturbance when the TEF case A control de- 
sign is implemented with analog components. The open- loop shock position response is 
included as a reference. Compared with the shock position feedback system shown in 
figure 15, low frequency attenuation is not quite so good. However, as frequency in- 
creases, the TEF system does better in the 6 to 20 hertz range. Since there is less lag 
between the p te signal and the disturbance than between y g and the disturbance, it is 
expected that this system might have an easier job of attenuating disturbances at the 
higher frequencies where phase lag is becoming a problem. Beyond 20 hertz, however, 
the system appears open loop. Again, this is because for case A little disturbance en- 
ergy exists in this region. 

Figure 19 compares the experimental and analytical frequency responses of the 
analog version of TEF case A control. Responses compare very well out to 20 hertz. 
After 20 hertz the comparison is not good. The probable cause is that beyond this fre- 
quency the analytical inlet model used for design and prediction was not an extremely 
close fit on amplitude. This can be seen by looking at figures of the curve fit informa- 
tion of figure 3 in appendix C. A more accurate fit would probably have produced closer 
agreement between experimental and analytical results. 

Figure 20 is a comparison of the closed- loop shock position frequency response of 
the analog and digital computer versions of the TEF case A control design. The two 
responses are almost identical. As stated earlier, the digital control algorithm was 
designed for and used a sampling time of 1 millisecond. This is quite adequate for dis- 
turbance frequencies up to 100 hertz; therefore, the close correlation with the continu- 
ous analog version as shown in figure 20 is as expected. 

Figure 21 presents the closed- loop frequency response of the analog TEF case B 
controller design. Also shown is the inlet open- loop frequency response. It should be 
remembered that case B has a certain amount of disturbance energy in the midfrequency 
range and not as much at the very low frequencies. As a result, the controller perform- 
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Figure 17. - Comparison of closed-loop experimental frequency responses of shock 
position to disturbance airflow using SPF case A analog and digital controls. 



Figure 18. - Comparison of experimental open- and closed-loop frequency responses 
of shock position to disturbance airflow using TEF case A analog control. 
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Figure 19. - Comparison of analytical and experimental closed-loop frequency 
responses of shock position to disturbance airflow using TEF case A analog 
controls. 



Frequency, Hz 

Figure 20. - Comparison of closed -loop experimental frequency responses of shock 
position to disturbance airflow using TEF case A analog and digital controls. 
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ance shown in figure 21 does not attenuate the low frequency disturbances as much as 
either of the case A controllers, but it does produce more attenuation in the midfrequen- 
cies out to about 20 hertz. This is shown by figure 22 which compares experimental 
analog versions of the TEF system for both the case B and case A designs. Case A 
provides very little attenuation after 3 to 4 hertz, whereas the case B design '’keeps 
working" out to 20 hertz. The slight magnification over the open loop shown in figure 21 
is difficult to explain except that control in this region is quite difficult because of the 
great deal of phase lag from the inlet at these frequencies. 

Figure 23 presents a comparison of the experimental response of the analog TEF 
case B design and its analytical counterpart. The prediction is quite good, especially 



Figure 23. - Comparison of closed-loop analytical and experimental frequency 
responses of shock position to disturbance airflow using TEF case B analog 
controls. 


in the midfrequency range. In the area of 80 to 100 hertz the analytical prediction shows 
better disturbance attenuation than the experimental. This is probably due to the limited 
order inlet model used in the analytical design. The additional phase shift of the actual 
inlet causes the experimental data to be somewhat degraded in this region. 

Summary of experimental results . - To summarize the experimental data, several 
observations can be made. In general, the three controller designs performed in agree- 
ment with their analytical predictions except for the shock feedback system SPF, where 
a difference in inlet model gain was found to exist. Each of the three controllers was 
implemented with an analog computer. Only two of the three were put into discrete form 
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and implemented with a digital computer. Where the two were compared to their analog 
counterparts, equivalency was quite good. In general, the experimental program point- 
ed out the problems of designing and then implementing linear optimal controller de- 
signs. 

Conclusions and Recommendations 

It has been demonstrated that inlet controllers which minimize the expected frequen- 
cy of unstarts can be designed and implemented. It was shown that controller charac- 
teristics depend strongly on the spectrum of the disturbance. Both shock position (SPF) 
and throat exit static pressure (TEF) were used as feedback variables. Because of the 
difference in noise levels on these signals, it was found that SPF was better for disturb- 
ances rich in low frequencies and TEF was better for higher frequency disturbances. 

Analytically predicted and experimental closed- loop frequency responses were found 
in general to be in close agreement. Experimental controllers were implemented with 
both an analog computer and a digital computer. Analog and digital results compared 
quite favorably. 

This attempt at applying linear stochastic optimal control theory to inlet control 
problems has, therefore, met with some success. The possible benefits which can be 
gained from using this relatively new theory seem to be great. However, this investi- 
gation exploited only a small fraction of the capabilities of the stochastic optimal control 
design approach. The remainder of this section contains recommendations as to areas 
that might warrant further investigation, both analytically and experimentally: 

(1) Consideration should be given to designing an optimal controller for both shock 
position and throat Mach number for both atmospheric and compressor face disturb- 
ances. This problem could use the multiloop capabilities of the optimal control design 
technique. 

(a) Multiple measurements could be used: throat exit and compressor face 
pressure, throat Mach number, cowl lip Mach number, etc. Experiments would 
be needed to determine more precisely various measurement channel noise levels 
and spectrum of the compressor face disturbance. 

(b) Controlled variables such as spike position (and possibly engine speed) 
as well as bypass door opening should be considered. 

(c) A quadratic performance index could be used involving mean- square 
values of shock position and throat Mach number plus penalties on bypass door 
opening, spike position, actuator slewing velocities, etc. 

(2) Sensitivity studies should be conducted on future inlet controller designs to de- 
termine the degradation in performance when inlet or noise parameters vary from those 
assumed in the design process. 
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(3) Controllers developed in the present study tended to be somewhat complex. 
Therefore, methods of developing simpler optimal or suboptimal controllers should be 
studied. Some approaches might be as follows: 

(a) Reduce the order of the open- loop inlet model and compare the results 
(on the more complete inlet model) using controllers based on the lower order 
model with those of the more complex model. 

(b) Investigate techniques of simplifying controllers which have been designed 
for the complete inlet model. 

(4) The approach used in this study for digital control is not unique; thus, alternate 
approaches to optimal digital computer control of inlets should be studied. The goal is 
to increase the sampling period while achieving performance comparable to continuous- 
type control. This could be done by: 

(a) Studying various ways of discretizing inlet and/or controller differen- 
tial equations. 

(b) Developing a method for directly determining the optimal discrete-time 
control law for a continuous-time system. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, November 10, 1972, 

501-24. 



APPENDIX A 


SYMBOLS 


A 

B 

C 


D 

d 

G BPD (s) 

G DUCT^ s ^ 

G INLEt( s ) 

g ns (s) 

G SHOCK^ s ^ 

H 



m 



system matrix, n x n 
control matrix, n x c 
output matrix, ixn 
dimension of u 

plant disturbance matrix, nxd 
dimension of w 

bypass door transfer function, (kg/sec )/V 
inlet duct transfer function, (N/cm )/(kg/sec) 
overall inlet transfer function, cm/(kg/sec) 
noise shaping transfer function, ND 

9 

inlet shock position transfer function, cm/ (N/cm ) 

measurement matrix, m x n 

identity matrix 

performance index 

equivalent quadratic index 

control gain matrix, c x n 

estimator gain matrix, n x m 

weighting factor 

dimension of y 

dimension of z 

dimension of x 

diagonalization transformation matrix 
2 

pressure, N/cm 

2 

throat exit static pressure, N/cm 
transformed state vector, n x m 
number of terms in truncated series 



s Laplace variable, sec" 

T sampling period, sec 

t time, sec 

u control vector, c x 1 

v measurement noise vector, m x 1 

v shock position measurement noise, cm 

y s 2 

v throat exit static measurement noise, N/cm 
p te 

Wj equivalence coefficient 

Wg equivalence coefficient 

w plant disturbance vector, d x 1 
w control airflow, kg/sec 

V 

compressor face disturbance airflow, kg/sec 
Wj total inlet airflow, kg/sec 

x state vector, n x 1 

y output vector, l x 1 

y_ shock position, cm 

Z arbitrary square matrix 

z measurement vector, m x 1 

a shock position tolerance, cm 



noise shaping transfer function parameters, rad/sec 

discrete estimator control input vector, n x c 
discrete measurement vector input to estimator, n x m 
discrete plant control input vector, n x c 
discrete disturbance vector, n x d 
expected frequency of unstarts, unstarts/sec 


<p ^ closed loop discrete state transition matrix, 2n x 2n 
<P e estimator discrete state transition matrix, n x n 
cp plant discrete state transition matrix, n X n 

Jr 
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a u RMS control effort, V 

mean-square control effort, 

2 

aj mean- square value of control airflow, (kg/sec) 
cfi mean- square value of disturbance airflow, (kg/sec 



2 

mean- square shock position, cm 

2 

mean- square shock velocity, (cm/sec) 
dummy variable 


r d transportation lag, sec 

Superscripts: 


differentiation with respect to time 


trial values 


optimal estimate of a vector 
T transpose 

Operators : 

PSD( ) power spectral density of 
( ) j normalized to 1 Hz 
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APPENDIX B 


APPARATUS AND PROCEDURE 
Inlet Description 

The inlet used for the investigation was an axisymmetric mixed compression type 
with 60 percent of the supersonic area contraction occurring internally at the design 
Mach number of 2. 5. A cutaway view of the NASA designed inlet is shown in figure 24. 
Specific characteristics of | the inlet as well as the tunnel test conditions are tabulated in 
table IV. Additional aerodynamic design details and steady- state performance charac- 
teristics of the inlet are given in references 21 and 22. The dynamic responses of vari- 
ous inlet internal pressures and of normal shock position to airflow disturbances are 
reported in reference 6. 

Shown in figure 24 are the inlet's translating centerbody and overboard bypass doors 
of which there are a total of six. Both the bypass doors and centerbody are hydraulically 
actuated and electronically controlled. Three of the symmetrically located bypass doors, 



Figure 24. - Cutaway view of inlet model. 
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TABLE IV. - DETAILED INLET SPECIFICATIONS 


AND TUNNEL TEST CONDITIONS 


Inlet 

Cowl lip diameter, cm 

47.3 

Capture area, cm 2 

1760 

Capture corrected airflow, kg/sec 

16.2 

Tunnel test conditions 

Mach number 

2.5 

Total temperature, K 

315 

Total pressure, N/cm 2 

8.95 

Specific heat ratio 

1.4 

Reynolds number (based on cowl- lip 

3.88XI0 6 

diameter) 


Angle of attack, deg 

0 

Inlet orifice termination description 

o 

Choke plate area, cm 

598 

Flow coefficient 

0.985 

Location of choke plate from cowl lip, 

146.5 

cm 



driven in parallel, were used to provide sinusoidal disturbances in diffuser exit correct- 
ed airflow. The remaining three bypass doors, also driven in parallel, were used as the 
manipulated variable of the various normal shock controllers. 


Inlet Instrumentation 

Figure 25 indicates the location of pressure taps connected to dynamic strain gage 
pressure transducers used in the investigation. The pressure transducers were close 
coupled to the pressure taps to enhance their response capabilities. Details of the 
location, response, and usage of the pressure sensors are documented in references 
8 and 16. 
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Station, cm from cowl lip 

Figure 25. - Details of inlet including instrumentation locations. 

Shock Sensing 

In this program eight throat static pressure signals were used as inputs to an elec- 
tronic normal shock position sensor. The logic required for this sensor was imple- 
mented on both a general purpose analog computer and on a digital computer. The de- 
sign details of this sensor are discussed in references 12 and 16. The output produced 
by either implementation was a stepwise continuous signal indicative of shock position. 
The various shock position controls tested used either the throat exit static pressure 
Pte or stepwise continuous shock position sensor as the feedback signal for control. 

Controller Implementation 

The inlet controllers designed by modern control techniques were implemented on 
both a general purpose analog computer and on a digital computer. 

Figure 26 is a photograph of the general purpose digital computer used for imple- 
mentation of both inlet and engine controls. The system, shown in block diagram form 
in figure 27, consists of four major units: 

1. A digital computer with 16 384 words of memory, a read-restore memory cycle 

time of 750 nanoseconds, and a word length of 16 bits 

2. A digital interface capable of converting both analog and frequency signals to 

computer compatible digital words and converting computer generated words to 
analog and logical outputs 
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Figure 26. - Digital control computer system. 
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Figure 27. - Schematic of digital computer setup. 





TABLE V. - DIGITAL CONTROL COMPUTER SYSTEM CAPABILITIES 


Digital computer 


Magnetic core memory size, words . . ........ . 16 384 

Word length, bits plus parity 16 

Memory cycle time, nsec ................. . . . . .............. 750 

Add time, fisec 1.5 

Multiply time, jusec 4.5 

Divide time, jjsec ....................................... 8.25 

Load time, (jsec . . . .............................. 1.5 

Indirect addressing , . ............................ . Infinite 

Indexing . ....................... Total memory 

Priority interrupts .... . . . . . . . 28 separate levels 

Index registers 2 

Interval timers ........... . . . 2 


Analog acquisition unit 


Overall sample rate (maximum), kHz ................ 20 

Resolution of digital data, bits . . . 12 (plus sign) 

Output code Two’s complement 

Number of channels 64 

Input range, V full scale ±10 

Conversion time, jusec . 38 

Total error with calibration, percent . 0. 073 


Analog output unit 

Total number of digital -to -analog conversion channels (DAC) 
Resolution 13 bit DAC (10 channels), bits (12+1) ..... 

Accuracy (13 bit) DAC, percent of full scale ....... 

Resolution 12 bit DAC '(16 channels), bits (11+1) ..... 

Accuracy (12 bit DAC), percent of full scale ....... 

Output voltage range, V full scale 

Slew rate, V/jtxsec . . 


Priority interrupt processor 

Number of channels 

Input voltage range, V .... . 

Computer switching 

Comparator hysteresis, mV . . 
Comparator output, V 


10 

..... ±10 

Trigger on rise or fall 
. . . Adjustable from 35 to 650 
7 








3. A signal processing unit which provides signal conditioning and monitoring capa- 

bility between the digital interface and the propulsion system to be controlled 

4. Programming peripherals consisting of a high-speed, paper-tape reader and 

punch and a teletype 

The capabilities of the system are given in table V and a comprehensive description 
is available in reference 17. 

All inlet pressure measurements were passed through signal conditioners and iso- 
lation amplifiers to provide high-level (-10 to +10 V) inputs to the digital interface 
equipment. This unit contains a random access multiplexer, a sample and hold ampli- 
fier, and a 13-bit digitizer. The complete digitizing process from channel sample com- 
mand to entry of the digitized measurement into computer memory requires 50 micro- 
seconds. This process is automated through the use of a block data transfer unit which 
ties up the main frame for only one memory cycle per word transferred. Completion 
of the sampling process is conveyed to the computer by a priority interrupt from the 
block data transfer unit. 

Digital commands are issued directly from the computer main frame to the 13-bit 
digital- to- analog converters. These outputs are passed through isolation amplifiers to 
provide ground isolation of the digital system and then to the servoamplifiers driving the 
manipulated variables. 


Test Procedures 

Both open- loop (no input to control bypass doors) and closed- loop frequency re- 
sponse tests were run. For all tests, the steady-state operating point of the normal 
shock was located near the middle of the eight throat static pressure taps. This was 
accomplished with an appropriate steady- state setting of the six bypass doors. An ap- 
propriate disturbance amplitude was determined from the open- loop response tests. For 
the open- loop tests, the three disturbance doors were oscillated sinusoidally at an am- 
plitude sufficient to move the normal shock over the eight throat static taps (fig. 25) at 
1 hertz. This was the disturbance amplitude used at all frequencies for both the open- 
and closed- loop testing. In addition to open- loop tests, closed- loop frequency response 
tests were run using both the throat exit static pressure p te and the stepwise continuous 
shock position sensor output as measured variables . These tests were intended to de- 
termine the capability of the feedback controllers to regulate inlet shock position in the 
presence of compressor face disturbances. 

For the frequency response tests, both magnitude and phase data for a few signifi- 
cant signals were determined online using a commercial frequency response analyzer in 
the control room. These signals, as well as many others, were recorded in analog form 
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on magnetic tape for reduction at a later time. All frequency response data were plotted 
in the form of Bode plots. The magnitude response data were normalized to the magni- 
tude at 1 hertz. 



APPENDIX C 


FREQUENCY RESPONSE CHARACTERISTICS OF THE UNCONTROLLED 
(OPEN-LOOP) 40/60 SUPERSONIC INLET 

A general physical description of the mixed compression (40/60) inlet used for the 
program discussed in this report is presented in appendix B (figs. 24 and 25). Frequen- 
cy response data describing the dynamic characteristics of this inlet were obtained in 
the test programs described in references 7 and 8. Figure 28 describes in block dia- 
gram form the particular inlet transfer functions which were obtained. Figure 28 indi- 
cates that an incremental airflow disturbance Wj occurring at the downstream or com- 
pressor face end of the inlet will propagate upstream, resulting in a pressure variation 
at the throat exit (p te ) pressure sensor location. It will also cause motion of the normal 

shock y 0 about its quiescent or desired steady- state location, 
s 



g inlet ( s) * ^ (s) 

Figure 28. - Block diagram of inlet characteristics to downstream airflow disturbance. 

The data presented in this appendix are used as a basis for determining transfer 
function relations G DUCT , Gj NLET , and G gH Q CK which approximately describe the 
inlet. These relations serve as a starting point with which to undertake the control de- 
sign in the CONTROLLER DESIGN AND ANALYTICAL PERFORMANCE section. 

Shock Position Measurement Model (G||\|lej(s)) 

The first configuration to be considered is the overall inlet response of shock posi- 
tion y„ to the disturbance w- . The data points in figure 29 show the experimental am- 
plitude and phase frequency response performance of the inlet shock position y g to the 
airflow disturbance w i at the compressor face station. The amplitude data have been 
normalized to the value at 1 hertz. To ensure the linearity required for future transfer 
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Figure 29. - Comparison of analytical and experimental frequency response performance of shock position to 
inlet airflow. 


function representations of the inlet frequency response data, the airflow disturbance Wj 
was of a small enough amplitude to encounter as few nonlinear effects as possible. 

An analytical representation of the experimental data presented in figure 29 was ob- 
tained by curve fitting the frequency response characteristics of an approximate transfer 
function model to the amplitude and phase data. Equation (Cl) is the result of the curve- 
fitting effort: 



Included in this equation is the steady-state gain relation for the inlet. The solid lines 
of figure 29 give the frequency response of equation (Cl), where the amplitude response 
has been normalized to the amplitude response at 1 hertz. The approximation of equa- 
tion (Cl) is quite good out to 100 hertz for both amplitude and phase. Beyond 100 hertz 
the amplitude response of the approximation falls off and deviates from the experimental 
data while the phase angle is still reasonably accurate. The approximation could, of 
course, be improved by the addition of more poles and/or zeros to the transfer function 
Gjnlet^ e Q. ua ti°n (Cl). It is felt, however, that improved higher frequency (100 to 
200 Hz) model accuracy at the expense of increasing the complexity of equation (Cl) was 
not warranted. 


Throat Exit Static Pressure Measurement Model (Gquq(s) and G^j^qq^s)) 


The second configuration to be considered is a model which involves two distinct 
frequency response relations. These are the relations of the shock position measure- 
ment to a variation in the throat exit static pressure and the variation of p^. e to a com- 
pressor face disturbance Wj. 

The data points of figure 30 show the experimental frequency response of the throat 
exit static pressure p tg to the disturbance w.. The coefficients in the transfer func- 
tion of equation (C2) were found by curve fitting the data of figure 30: 


J DUCT 


1. 52 - 

(s).fss(s). m 


+ 1 


s A -1.5X10" 3 s 
+ lie 


\500 


w i 


+ 1 


s V , 2(0. 3)s 


N/cnot 

kg/sec 


(C2) 


^80 


1365, 


365 


+ 1 
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The frequency response performance of equation (C2) is given by the solid lines of 
figure 30. The amplitude response fits well to about 80 hertz at which point it drops off 
and does not duplicate the resonances of the experimental inlet data. The phase angle 
response in figure 30(b) is quite good to 140 hertz. The approximation therefore is felt 
sufficiently accurate for the purpose of controls design. 

The data points in figure 3 1 show the experimental frequency response of shock 
position y g to the pressure p te - Since all frequency response information was deter- 
mined by an airflow disturbance Wj at the compressor face, the data in figure 31 were 
obtained by finding the difference between the experimental data in figures 29 and 30. 
The transfer function approximation which was curve fit to the data in figure 31 is given 
by 


g shock 


(s) — — — (s) = 
Pte 


10. 68e 


•2.5xl0“ 3 s 



cm 

N/cm 3 


(C3) 


The frequency response performance of equation (C3) is given by the solid curves 
in figure 31. 
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Figure 31. - Comparison of analytical 
throat exit pressure. 
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APPENDIX D 


DISCRETE CONTROLLER FORMULATION 

The inlet discussed in this report is being controlled by an optimal feedback con- 
troller described by continuous linear time-invariant differential equations. When ex- 
pressed in matrix form, these are equations (16) and (17). These equations, when they 
have been slightly rearranged, are 

£(t) = (A - K e H)x(t) + Bu(t) + K e z(t) (Dl) 

u(t) = -K c x(t) (D2) 

Note that z(t) is the controller input and u(t) the controller output. As discussed in the 
main text, the optimal controller was implemented with a digital computer as well as an 
analog computer. The digital version is shown in figure 14. In order to obtain a digital 
controller algorithm, a discrete-time approximation must be obtained for the continuous- 
time equations (Dl) and (D2). The method used in this report is discussed in the follow- 
ing sections. It should be pointed out that the method used for obtaining the digital con- 
trol algorithm is not unique, and the technique used is only one of a number of possible 
approaches. 


Discrete Control Algorithm 


The solution to the vector- matrix differential equation (Dl), at time t^ + p given the 
state x at time t^ is 


* (t k+l } = ^e (t k+l " + / tk+1 ^e (t k+l “ t)Bu(t) dT + / tk+1 ^e (t k+l " T)K e z(r) dr 

*k *k 


(D3a) 


where 

*e<W ' V = ex P[< A ' K e H ><*k + l ' *k>l 

’’e^+l ' T) = exp [ (A ' K e H)(t k + l - T >] 
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Let = kT and = kT + T be successive sampling instants, separated by sample 
period T, If it is assumed that u and z are constant over kT s t < kT + T, z and 
u can be moved out from under the integral signs in equation (D3a) and an expression 
can be obtained for x(kT + T) in terms of x(kT), u(kT), and z(kT). Since the digital 
computer produces a stepwise continuous signal u, u is constant during a sample time 
(u(r) = u(kT), kT s T<kT + T), However, z is not; thus, we must assume that T is 
small enough so that z is approximately constant during the interval: 

z(r) = z(kT) kT < t < kT + T 


Making these substitutions, equation (D3a) becomes 

f rkT+T I 

x(kT + T) = ^ e (T)x(kT) + jjf ^ e gk + 1)T - t]B drV u(kT) 


/•kT+T | 

J 9 e [(k + 1)T - r]K e drl z(kT) (D3b) 
kT " J 


But since the two integrals in equation (D3b) are independent of k, we can evaluate them 
for k = 0: 


x(kT + T) = <? e (T)x(kT) + 


x(kT + T) = (p (T)x(kT) + 


rT 


/-T I 

1 f p (T - r)BdT 
-0 e 

u(kT) + 

J g ^ e (T-T)K e dr 


r T 1 


" /- T 

1 cp A t) B dr 

L/o J 

u(kT) + 

J 0 '?e (T > K e dT _ 


z(kT) 

(D3c) 

z(kT) (D4) 


Since 

(A-K H)t 

cp e (r)=e 


then 


d <p (A-K H)t 

£ = (A - K 0 H)e e = (A - K e H)<P e (T) 

dT 

Therefore 
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<pJj) dr = (A - K.H)’ 1 d<p 


(D5) 


and 


r T t 

T -1 -1 

[ (A-K H)T 1 

1 <pJt) dr = I 
'0 e J ( 

^ (A - K e H) 1 dcp e = (A - K e H) 1 

e e - I 


(D6) ' 


Substituting equation (D6) into equation (D4) yields 


r -i 

f (A-K H)T ]1 

j(A - K e H) 1 

P "JJ 


kBu(kT) 


+ UA- K e H) 


-1 


(A-K H)T 


- I 


►K e z(kT) 


Therefore, the discrete controller computer algorithm is defined by equations (D7) and 
(D8): 


x(kT + T) = <pj T)x(kT) + r ce u(kT) + r m z(kT) 


m 


(D7) 


where 


<P e ( T) = e 


(A-K H)T 


r ce=< A - K e H >' 


(A - K H)T 


- I 


B 


r m = ( A " K e H > 


-1 


(A-K H)T 


- I 


K 


and 


u(kT) = -K x(kT) (D8) 

The matrices <p (T), T _ . , and F „ must be numerically determined for the appro- 

G CG XU 

priate sampling time T. Acceptable sampling times are those which result in accept- 
able stability of the complete system operating with the sampled-data controller. 
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Stability Considerations 


To determine the stability of the closed- loop sampled- data system, the continuous 
inlet plant described by equations (3) and (5) (repeated here as eqs. (D9a) and (D9b)) can 
be put in discrete form : 


x(t) = Ax(t) + Bu(t) + Dw(t) 

(D9a) 

z(t) = Hx(t) + v(t) 

(D9b) 

x(kT + T) = (p v ( T)x(kT) + T cp u(kT) + T d w(kT) 

(DlOa) 

z(kT) = Hx(kT) + v(kT) 

(DlOb) 


where 


<P p (T) = e AT 

r cp =A" 1 (e AT -I)B 
r d = A" 1 (e AT - I)D 


Here again, an approximation is made that w(t) and v(t) are assumed to be constant 
over the sampling interval. Combining equations (D7), (D8), (DlOa), and (DlOb) results 
in the following expanded matrix equation, which represents the closed- loop sampled- 
data system: 


~x(kT + T)~ 

xfcYVrT 


>e (T) - r ce K c ! r m H ' 


"x(kT)~ 


"r ! o" 

m I 

, . j, 


"v(kT)~ 

- r cp K c ^p(T) 


_x(kf)_ 


_° ! r d 


w(kT) 


(Dll) 


Let 


.> e (T) ‘ r ce K c r m H 1 
<p (T) = - e _c i _m_ 

CL 1 r„„K„ ]<p{ T) 


cp c 


i 


(D12) 


The eigenvalues of <^ ( ~, L (T) determine the stability of the closed- loop system. A 
sampling time T which results in an eigenvalue with a magnitude greater than one 
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(outside unit circle in the complex plane) is unacceptable. This criteria was used in 
determining the largest acceptable sampling period. 


Matrix Exponential Expansion 

(A-K H)T AT 

Matrix exponentials such as e and e can be determined by using a 

matrix form for the series expansion of e raised to some power. This technique was 
used in arriving at an acceptable method for determining the matrix coefficients of equa- 
tion (D7). A brief description of the procedure using a square matrix Z is 


e Z T =I+ZT+ zV + zV + . 


2 ! 


3! 


(D13) 


and 


Z-^e 21, 


I) = Z" 1 


4t + 



(D14) 


If the series is truncated after r terms, equation (D14) becomes 



r 


> 


T + ^T 

T + . . . + ZT (T + ZT V 

> 

2 

3 

r-A r / 




\ / J 

J 


(D15) 


and equation (D13) becomes 


e ZT = I + Z< 


T + 


ZT 


T + 



(D16) 


The computer subroutine STM whose listing appears at the end of this appendix actually 
implements equations (D15) and (D16). Enough terms of the expansion are used to en- 
sure that the matrix elements have converged with sufficient accuracy. 


Numerical Considerations 

The numerical values of the elements of the matrices <^(T), r , and r for a 

C Cc XI* 
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stable sampling time T were found to have a wide spread. The computer on which the 
control law of equations (D7) and (D8) was programmed used fixed- point machine lan- 
guage. Thus, the numerical spread of the matrix elements created significant scaling 
problems in programming the controller. Also, the vector-matrix multiplications re- 
quired to implement equation (D7) (especially <? e (T)x(kT)) took too much computer time. 
This was because <P e (T) had few nonzero elements. To alleviate some of these prob- 
lems, steps were taken to (1) condition the numerical elements to reduce scaling prob- 
lems, (2) reduce the number of elements in the <p 0 matrix, and (3) provide a check on 
the final results. 

It was found that a convenient way to accomplish steps (1) and (2) was to use a block 
diagonal transformation on <p (see ref. 23), This brought the numerical values of <p a 
closer together and eliminated many of the multiplications required in executing the 
computer control law. The block diagonal transformation is now outlined. Let 


Pq(kT) = x(kT) 

(D17) 

Pq(kT + T) = x(kT + T) 

(D18) 


Substituting equations (D17) and (D18) into equations (D7) and (D8) yields 
Pq(kT + T) = <p e (T)Pq(kT) + r ce u(kT) + r m z(kT) 


or 

q(kT + T) = P' 1 ^ (T)Pq(kT) + P _ 1 r n u(kT) + P -1 r z(kT) (D19) 


and 


u(kT) = -K c Pq(kT) (D20) 

The matrix P is a transformation matrix whose columns are the eigenvectors of 
cp e { T) if all the eigenvalues are real. If there exists a complex conjugate pair of eigen- 
values, the column of P which would correspond to the first half of the eigenvalue pair 
is set equal to the (vector) sum of the pair of eigenvectors. The column of P which 
would correspond to the second half of the eigenvalue pair is set equal to the difference 
of the eigenvector pair. The resulting block diagonal matrix, P~^cp (T)P has the real 

V 

eigenvalues of <^ e (T) on the diagonal. When complex conjugate pairs are present, the 

real parts lie on the diagonal with the imaginary parts on the off- diagonals. Where « (T) 

2 - 1 e 
has n nonzero elements, P <p e (T)P has n, if all eigenvalues are real, and 3n - 2 
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if all eigenvalues are complex. This represents a considerable reduction in computer 
operations required in implementing equation (Dip) as compared to equation (D7). 

Once the transformation matrix P is determined, the numerical values required by 
the controller can be determined. The transformation results in a transformed set of 
estimated states q. The control output is determined by the transformed control 
weighting -K C P on these new estimated states. In determining the control output, each 
of the transformed state estimates must be calculated. Therefore, the programmer 
must know what will be extremes or worst case values of these states when operating 
within the closed- loop experimental system so that he can properly scale the states. To 
assist the programmer in this area, the system was analyzed analytically to determine 
worst case magnitudes of the estimated states. This was done by subjecting the closed- 
loop discrete system to inputs equivalent to the worst case compressor face airflow dis- 
turbance anticipated for the experimental program. For the transformed controller, 
the closed- loop system, as defined by equations (D19) and (D20) together with plant 
equations (DlOa) and (DlOb), is 


q(kT + T) 

] 

x(kT + T) 



<? e (T)p-p- l r ce K c P jp ' 1 

” ^ cp K ? j 


r h 

m 




p _1 r ' o 

m I 


V(kT)" 

w(kT)_ 


(D21) 


The computer subroutine used to accomplish evaluation of the transient performance 
of the closed- loop sampled- data inlet system is included in this appendix. 


Supporting Computer Programs 

In this section is a brief description of the computer routines used in a large central 
batch processing facility to arrive at an acceptable discrete controller capable of being 
expeditiously and reliably programmed on the digital control computer system . 

Shown in figure 32 is a flow chart of the computer program and its associated sub- 
routines. Once the appropriate constants representing the continuous inlet and its 
selected estimator/ controller are read into the computer, the programs use subroutine 
STM and a preselected value of sampling time to determine the discrete equivalent of 
the control and inlet plant. Eigenvalues of the <p ^ matrix are then determined for 
the particular T used. If the eigenvalues are all within the unit circle, then subroutine 
DIAG is used to effect a coordinate transformation of the state estimator. Finally, the 
complete closed- loop sampled -data system is exercised through a transient. These 
transient data show the sizes of the various computer estimated states and feedback con- 
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Figure 32. - Flow chart for digital computer program (DIGCON). 








trol input and thus are an aid in machine language programming the computer controller. 

Listings of the computer programs used in determining the discrete control laws 
are included in this appendix. Table VI lists the significant equation variables with their 
corresponding FORTRAN designations. Other information concerning the attached pro- 
gram listing is contained in various comment statements appropriately located through- 
out the program. 


TABLE VI. - CROSS REFERENCE LIST 
FOR DISCRETE CONTROLLER 
VARIABLES 


Equation variables 

FORTRAN variables 


PHE 

r ce 

GAMU 

r m 

GAMZ 


PHP 

p cp 
r d 

GAMB 

GAMD 

K c P 

KCP 

p" 1 .r m 

GAMZP 

p _1 iv 

ce 

GAMUP 

p _ 1 r„ k p 

ce c 

GUKC 

p _ 1 r m H 

m 

GZH 

r G P K C P 

GBKC 




ilBFTC DIGC.CN LIST 
C 

C **#*#***#******$********>!:#*****##*#*****#******##******#*##********#** 

c 

C f) I GCON DETERMINES A DISCRETE FORM OF A FEEDBACK CONTROL 

C L Am SUITABLE FuR IMPLEMENTATION WITH A DIGITAL COMPUTER. 

C CONTROL WEIGHS ESTIMATES OF THE SYSTEM STATES OBTAINED WITH 

C A KALMAN FILTtk. PREDICTION CF TRANSIENT RESPONSE PERFORMANCE 

C OF CLOSED LOUP SYSTEM TO A STtP DISTURBANCE IS ALSO INCLUDED. 

C 

C THE FOLLOWING SL8RGUT INE S AkE PART CF TFE IBM SSP 

C GMPRD - MULTIPLIES TWU GENERAL MATRICES TO FORM 

C A RESULTANT GENERAL MATRIX. 

C FACTR - PROVIDES FACTORIZATION OF THE INPUT MATRIX 

c into a product of a lower triangular matrix 

C AND AN UPPER TRIANGULAR MATRIX. 

C HSBG — REDUCES A REAL MATRIX TO ALMOST TRIANGULAR 

C (HtSSENbURG)FQkM. 

C DMINV - INVERTS A DOUBLE PKECISILN MATRIX YIELDING 

C A DOUBLE PRECISION RESULT. 

C 

C OTHER SUBROUTINES CALlEU BY UIGCON ARE E I GOR » STM. AND DIAG. 

C LISTINGS FOR THESE ARE INCLUDED. 

C. 

C Jit******#*# ##*##########*##*######*#*##*######*####* 

f 

C DIMENSION OF VARIABLES 

C A ( N , N ) 

C HlN.C) 

C H ( M . N ) 

C KCIC.Nl 

C K F I N . M i 

C OtN.C) 

C CC(M.N) 

C SING(N.N) 

C 1 . BkCIN.NI 

C GUKC(N.N) 

C FXTl(N.N) 

C FCN.N) 

C PHP(N.N) 

C GZHIN.N) 

C GAMAT(N.N) 

C GAMBIN.C) 

C GAMCIN.C) 

C KCP(C.N) 

C XS(M« 1 ) 

C P56IM.1) 

C PFfllGIFN.PN) 

0 PHFSfi ( 2N . 2N ) 

C X B I G ( 2 N » T > 

C XNBIGI2N. I > 

f. / B I G ( ? . 1 ) 

c fxibk;i2n. 1 > 

C EX2rtlG(2N.l) 

C GAMBlG (2N.2) 

f RR2(2N) 

C RI2I2N) 

C MAG2I2N) 

C PFE(N.N) 


62 



n ft n n ^ n 'h n n n n ft n ^ rj n 'ft h n n n n n n 'n h 


c gamu(n.C) 

i: gamz(n.c > 

C TC(N,N> 

C PMO(N.N) 

C GAMUP(N.C) 

C GAM/PU.C) 

C. A A A ( N * N ) 

C FIG(W.N,2» 

f. CPR(N) 

C CPI IN) 

C RR(N> 

C R I ( N ) 

(. CB(N') 

C I ( N ) 

AH(N»N > 

i.wvn< n> 

MfcVD(N) 

PPHlD(N.N) 

PHIP(N.N) 

I T ( N • N ) 

EXT lb (M »N ) 

PFRN1N) 

1PFRMM 
l PER ( N ) 

DC(N.N) 

Fn(M.N) 

Y(N,1 ) 

S m N , 1 > 

FX(N. 1> 

U S 1 N » 1) 

JUOQdV.l) 

EXT2(2N,2N) 

PERPM2N) 

IPER2M2N) 

IPER2( 2N> 

P2S12N.2N) 

P212N.2N) 

02(2 N. 1) 

/ 2 ( 2N « 1) 

0S1< 2N.1 ) 

QP2(2N,1) 

yo2(2N*i) 

CLMMON / RLCA / A(IO.IO). B(10,i), HU. 10), KC(l.iO) * KE(lO.l), 

1 0(10,1). CCU, 10) 

DIMENSION SINGUO.iO), GBKU10.10), GUKCU0.10), 

1 EXT 1 ( 10.10). filO.10). P HP (10,10), GZHUO.iO), GAMATU0.10), 

2 GAMB ( 10, 1 ) , GAMDUO.l), KLP ( 1. 10) , XS(l.l), P56U.1), 

3 PHB I G ( 20 , 20 ) , PHESB ( 20 , 2 0) . XBIG120.1). XNBlGUO.l) . ZBIG(2,i), 

9 EX IB IG( 20, 1 ) , EX2bIG(20, 1) , GAMB I G ( 20, 2> » RR2(20), RI2(20), 

5 MAG2120), PHt ( 10.10). GAMOllO.l), GAMZUO.i), TOUO.iO), 

6 PHIC(IO.IO). GAMUP (10,1), GAMZPUO.l), AAAUO.IO). E IG( 10, 10 , 2) , 

) CPR(IO), CP I (10) . HR ( 10 ) , RKLO), CRUC), CIUO), AHU0.10), 

8 LwVDUO), MwVU(lO). PPHIU( 1U.10) , PH1PU0.10), TTU0.10). 

9 EXT 1 6 (10 • 10 ) , PERN(IO), IPERN(IO), IPER(IO), 00(10,10), 

A EO(IO.IO), YUO.l), SOUO.l) , EX (10, 1 ) , OSUO.l), 000(10,1), 

B ExT 2(20,20), PER2N(20), IPER2M20), IPER2(20), P2S(20,20), 

C P2(20,20>. 02(20,1) , 22(20,1), 0S1 (20. 1) , 0P2(20,1), YU2(20, 1) 

D01JBL E PRECISION A, B, H, KC , KE. EXT1. TD, F, DT , PHE, GAMU, 

1 GAMZ , PHP, bAMB. GAMC, 0. CL, PHIO, GAMUP, GAMZP, KCP, TT, PPHID, 

2 PHlP, EXIBIg, tX2b l G , X B 1 G , ZBiG, XnBIG, GAMBIG, P3b, XS, GAMAT, 
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G 

C 


f. 


C 

c 


G 

G 

G 


C 

C 

c 


c 


c 

c 


c 

c 


3 PHBIG, GZH. G6KC, GUKG 
REAL PAG* 

INTEGER C 

CfftMON / FORM / VFMIUO.bi. FMT (6) 


DATA ((VFMTCt.JJ, 

J 

= 1 9 O ) t 

I * i 9 10 ) 




1 

/ 6 H 

* OH 

* 

OH 

* 6H 

6 H 

i iPi. 

6HE12.4) , 

2 

6 H 

* 6h 

9 

6 H 

, 6H 

6 H 

{ iP2 , 

6 Hi: 12*41 1 

A 

6 H 

t OH 

* 

oH 

» oH 

oH 

1 1 P 3 . 

6HE12.4) » 

H 

6 H 

* 6n 

• 

6 H 

♦ oH 

6 H 

1 1P4, 

OHE12.4) t 

b 

6 H 

v OH 

• 

OH 

« 6H 

6 H 

(IP5, 

0HE12.4) t 

6 

6 H 

* 6H 

t 

OH 

, OH 

6 H 

( iP6, 

6HE12.4) , 

7 

6 H 

f Oh 

* 

6 H 

t 6H 

oH 

UP?, 

OHE12.4) , 

8 

6 H 

* OH 

t 

6 H 

♦ 6H 

6 H 

(IPS* 

6 HE i 2 * 4 ) f 

9 

6 H 

* 6H 

t 

6 H 

• 6H 

6 H 

i iP9, 

6 HE 1^*4) 9 

A 

OH 

« OH 

* 

OH 

« oH 

oH 

( 1P10, 

6 HE 12*4 ) 


IOP IS THE PRINT VARIABLE 

IF IOP = O. THERE RILL BE STANDARD OUTPUT 
IF IOP - 1. THERE WILL BE EXTENDED OUTPUT 

IOP = 0 
WRITE ( 6 . 35) 

H?. = ?G 
N = 10 

M = 1 
G = 1 
DT = .2D0 
KMAXE - 50 
K.MAXP = 50 

no 4 i=i » N 
on 3 k = i » c 
K.CPIX, II = 0.000 
OAMU(I.K) = C .000 
GAMZ(l.K) = 0.000 
GARB! I .K J = U.ODO 

3 G AMR I l «K 1 = O.UCO 
DO 4 J = l.N 
GUKC ( I » J ) = U.ODO 
GZhII.O) = O.ODO 
GRKCI I.J) = 0.000 

4 FXT 1 ( I.JI = 0.000 

FORM IKE * H) 


DO 5 K = l.N 
00 5 1 = l.N 

00 5 J = l.M 

5 FXU(I.K) = KEt I.J) * h(J.K) ♦ EXTIII.K) 

FORM (A - KE * F) ESTIMATOR MATRIX 

00 6 I = l.N 
00 6 J = l.N 

6 F C I.J* = A ( I.JI - E XT l ( I.J) 

IF imp .EU. 0 ) GO TO 329 


PRINT ESTIMATOR MATRIX 
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r> n r> n r> r> r> r> r> n n r n onn 


C 

WRITE (6.435) 

DC 501 J = 1*6 

501 FMT(J) ~ VFMT(N,J) 

WRITE (6.FMT) ((F(I.J), J * l.N), I * 1,N) 

329 00 330 I - l.N 

on 330 J - l.N 

330 S IMG < I . J ) = F ( I . J ) 

CALL FSBG { N. SING. N ) 

OfcTERMINE EIGENVALUES OF ESTIMATOR ( A - KE * H) 

IF ( IOP «EQ* 0) GO TO 3 26 
WRITE (6.475) 

328 CALL EIGQR (SING. N. RR. Rl* IOP) 

FLRM DISCRETE STATE TRANSITION MATRIX F OR ESTIMATOR ( A ~ K£ * H) 

CALL STM (F. PHE. GAMAT, N» KM AXE , OT, IOP) 

FORM AND PRINT THE INPUT VECTORS FOR THE ESTIMATOR 

00 331 K = i.C 

DO 331 I * l.N 

on 331 J - l.N 

GAMU(i.K) = GAMAT (I * J ) * B(J.K) * G AMU ( I * K ) 

331 GAMZ(I.K) * GAMAT (I .3) * KE(J.K) + GAMZ(I.K) 
w R I T E ( 6 , 35) 

WRITE (6 * 440 ) 

DC 502 J = 1.6 

50 2 EMT(J) = VFMT(C.J) 

WRITE (6.FMT) ((GAMO(I.J). J = 1,0* I * l.N) 

WRITE (6.45) 

WRITE (6.FMT) UGAMZ(I.J). J = i.C), I = l.N) 

FCRM DISCRETE STM FOR PLANT (A) 

CALL STM (A. PHP « GAMAT, N, KMAXP , DT , IOP) 

FCRM AND PRINT THE INPUT VECTORS FUR THE PLANT 

DO 332 K = i.C 

DC 332 I — l.N 

DO 332 J - l.N 

GAM 6 ( I * K ) - GAMAT ( I . J ) * BU.K) + GAM8(I,K) 

332 GAMO(l.K) = GAMAT ( I * J ) * D(J,K) + GAMOM.K) 

WRITE (6.35) 

WRITE (6,445) 

00 503 J = 1,6 

503 FMT(J) = VFMT(C.J) 

wRITE (6.FMTI ( ( GAM8 ( I * J) • J = 1,0, i = 1,N) 

WRITE (6.45) 

WRITE ( 6 , FMT ) ( ( GAMOCI ♦ J) ♦ J 1,0, I = l.N) 

DIAGONALIZE TEE DISCRETE STM FUR THE ESTIMATOR (A - KE * H) 

C 

CALL DIAG (PHE, G AMU , GAMZ , TD, PHIL, GAMUP, GAMZP, N, C, N2, AAA, 

1 EIG, CPR, CPI, CR, Cl, AH, LwVD, MwVD, RR, RI, PPHID, PH1P, TT, 

2 EXT 16 , PERN, IPtRN, IPER. DD, ED, Y, E XT 2 • SD, EX, US, UUQ, 

3 PEB2N « I P ER2, IPER2N, P2S, P2, U2, Z2, USi, UP2, YD2 , IOP) 
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C PRINT THE TRANSFORMED INPUT VECTORS FOR THE ESTIMATOR 

C 

WRITE 16.35) 

WRITE (6.450) 

DO 504 J * 1,6 
504 FMT(J) = VFMT(C.J) 

WRITE (6.FMT) 1 ( GAMUP( I, J) * J * 1,0, I * l.N) 

WRITE (6.45) 

WRITE (6.FMT) ( ( GAMZP ( I , J) , J = 1,0, I = 1,N) 

C 

C FORM AND PRINT THE TRANSFORMED CONTROL GAINS 

C 

DO 319 K * 1, N 
DO 319 I * 1,C 
DO 319 J = l.N 

319 KCPI I ,K) » KC(I,J) * TD( J »K) + KCP( I,K) 

WRITE (6,455) 

WRITE (6.311) 

DO 505 J * 1,6 
5(S5 FMT(J) as VFMT(N.J) 

WRITE (6.FMT) ( ( KCP ( J, I ) , I * 1,N), J = 1,0 

C 

C FCRM AND PRINT THE CLOSED LOOP DISCRETE STM 

C 

DO 405 J = l.N 
DO 405 I * l.N 
DO 405 K * 1,C 

GUKCI I , J ) * GAMUP(I.K) * KCP(K.J) + GUKC(I.J) 

GBKC( I.J) * GAMB(I,K) * KCP(K,J) + GBKC(I,J) 

405 GZH(I.J) = GAMZP ( I , K ) * H(K,J) + GZH(I,J) 

DO 410 I « l.N 
DO 410 J = l.N 

PHBIG(I.J) = PHlD(I.J) - GUKC(I.J) 

K a 1 + N 

PHBIG(K.J) a - GBKC(I.J) 

L a J + N 

PHBIG (1,0 a GZH( I.J) 

410 PHBIG(K.L) = PHP (I.J) 

IF ( IOP .EQ. 0) GO TO 411 
WRITE (6,35) 

WRITE (6.470) 

DO 506 J a 1,6 
506 FMT(J) = VFMT(N.J) 

WRITE (6.FMT) ( ( PHBIG! I.J), J = l.N), I a 1.N2I 

K = N + 1 

WRITE (6.FMT) ((PHBIG(I.J), J = K,N2), I = l,N2) 

411 DO 705 I = l.N 2 
DO 705 J = 1,N2 

705 PHESBI I.J) = PHBIG! I.J) 

CALL HSBG (N2, PHESB, N2) 

C 

C FORM AND PRINT EIGENVALUES OF CLOSED LOOP DISCRETE STM 

C 

CALL E IGOR (PHESB. N2, RR2, RI2, 0) 

WRITE (6,35) 

WRITE (6.460) 

WRITE (6,1002) 

00 415 I = 1.N2 

MAG2 ( I ) a SORT ( RR2( I ) * RR2(I) + RI2(I) * RI2ID) 
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r> n r; o n n r= r> n n n n r> n r> 


wKlTF (o. 10031 kR2 ( I > » RI2(I>. MAG21 1 ) 

415 f.nMINUE 

CHFCK STABILITY (if CLOSED LOOP SYSTEM 
EXIT IE UNSTABLE 

00 321 l = 1.N2 
Ir (MAG2<i) .GT. 1.0) GO TO 13 
321 CONTINUE 

DO 12 L L. = 1 . C 
00 322 I = 1,M 
P 36 ( I . 1 ) = O.UDO 
32 2 X S ( I * J ) = 0.000 
on 1 1 = 1 ,N2 

EX1HIGU.11 = U.ODO 
FX2B I G ( 1.1) = U.ODO 
XBIGI (.1) = C.OUO 
XNHIG(I.l) = U.ODO 
no 7 J = 1,2 

7 GAME l G ( I , J ) = O.OCO 

THIS IS A UNIT STEP 

/BIG (1,1) = O.ODU 

/HI G (2.1) = l.ODJ 

DO 8 I = i.N 

GAME I C- ( 1 « 1 ) = GAM/P ( I »LL ) 

K. = I + N 

8 GAME I G (K . 2 ) = GAMD(I.LL) 

COMPUTE OLUSEO LOOP STEP RESPONSE 

WRITE (6.35) 

00 12 K = 1,300 

on y i = i,n 2 
on ^20 J = 1 , N2 

420 FX18IGU.1) = PHBIG(l.J) * XBIG(J.l) + EXlBIb(I.l) 
00 9 J = 1,2 

S F X2H l G ( l. 1) = UAMB I G ( 1 * J ) * ZBIG( J, 1) + EX2BI G( 1,1) 
DO 10 I — i,N2 

X N B I G ( I . 1 ) = EXlttlG( l .l) + EX2b I G ( I , 1 ) 

10 XBIG(l.l) = XNE IG ( I , 1 ) 

PRINT ESTIMATOR STATES 

wR I T F (6,465) 

00 507 J - 1.6 
307 F«1(J) = VEMT ( N , J ) 

WRITE (6.FMT) (XNBIG(I.l). I = i.N) 

00 425 II = l.M 
DO 42 5 J * I.N 
JJ = J + N 

P 56 ( 11,1) = H(II.J) * XNBIG(JJ.I) + P56 (II, 1 ) 

425 XS(U.l) = CC(II.J) # XNE I G( J J , 1 ) + XS(il.l) 

PRINT PLANT UUTPUT RESPONStS 

WRITE (6,4 30 > 

WRITE (6,318) (P 36 (1,1), XS(l.i), I = l.M) 

DO 324 I = l.M 




P ifa ( I . 1 ) 

= 0 

.000 




X S ( t , 1 ) 

- 0 • 

ODo 




on i i i 

= 1 

. NR 




FX 1BIGU 

* 1 ) 

= O.UUO 




fxpbioi i 

♦ 1 ) 

= u.uOu 



1 1 

XNB I G ( i. 

1) 

0.000 



1^ 

( . L NT INUF 





3 5 

FO RMAT 

( 1 HI ) 



45 

FORMAT 

( 1 HO > 



1 JO? 

FORMAT 

( 1* 

/ LOX , 10H REAL PART, 10X, 11H 

I HAG* 

PART, iOXf 


1 10H MAGNITUDE) 



1003 

FORMAT 

( 8X . 

E14.7. 6X. E 14. ft o X, E14.7* 



J 1 l 

FORMAT 

{ bOX 

, SHKCP) 



318 

FORMAT 

( IX, 

IP2E20.O) 



-*30 

FORMAT 

( 1<*X 

, jHP 56, lax, 2HXS) 



435 

FORMAT 

( IX, 

IoHESTIMATUR MATRIX) 



440 

FORMAT 

( IX, 

SlHlNPUT VECTORS FUR THE ESTIMATOR) 


44 5 

F ( S k M A T 

< IX, 

27HINPUT VECTORS HOR ThE PLANT) 


*50 

FORMAT 

( lx. 

48hT RANSFORMEO INPUT VECTORS 

FOR THE 

ESTIMATOR) 

*5 5 

format 

l IX, 

28hTRANSFURME0 CONTROL GAINS) 



40 0 

format 

( IX , 

30Hfc IGENVALUtS OF CLOSED LCOP 

DISCRETE STM) 

46 5 

format 

tlx. 

22 HE ST I MATuR STATE VECTOR) 



*7 0 

FORMAT 

( IX, 

24FCLUSEU LO()P DISCRETE STM) 



475 

FORMAT 

(lx. 

24Ht IGENVALUtS OF ESTIMATOR) 



13 

CENT I NUF 

STOP 

FNO 
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* 1 FF TC STM LIST 

SUHKOliT l NE STM IF, PHI, GAMAT. N , KM AX « DT. IGP) 

r 

C **^ *#*#**#*#* **** $##****£** ** *#***£************* ****** 

C 

C STM COMPUTES THE STATE TRANSITION MATRIX PHI = EXP( F * OT ) 

C uSIUu A NESTED POLYNOMIAL b VALUAT I UN UF Thfc SEP-IEo REPRESENTATION. 

C 

C PHI IS THE OUTPUT S.T.M., GAMAT IS THE OUTPUT MATRIX FOR CUM- 

C PUTIN G THE VECTORS GAMMA. F IS THE INPUT MATRIX, N IS THE SIZE 

C UF THF MATRIX, kMAX IS THE NUMbER OF TERMS IN THE S.T.M. SERIES, 

C AND CT IS THE SAMPLING TIME. 

C 

C STM CALLS NO SUBROUTINES. 

C 

C ************* *** ******* **************** **** ** **** * * ** * ** * * *** * * ** * *** * 

c 

C DIMENSION UF INPUT VARIABLES 

C F { N » N 1 

C 0 IMFNS ION CF OUTPUT V Ah IABLES 

C PHI(NiN) 

t GAMAT ( N* N > 

DOUBLE PRECISION F, PHI, GAMAT , DT, X, Y 
DIMENSION FIN.N), PHl(N.N), GAMAT < N, N) 

COMMON / FORM / VFMT (10,0 , FMT ( 6) 

C 

c 

c 

C NN IS THF UNDERFLOW COUNTER 

C 

NN = 0 
C 

C INITIAL GAMAT I S I * CT 

C 

no ioo i = i , n 

DO 100 J = 1 * N 

gamati i, j > = o.odu 

IF (I .ED. J» GAMAT ( I, J) = DT 
100 CONTINUE 

KMAX1 = KMAX - 1 
C 

f. SERIFS CALCULATION OF GAMAT 

C 

DC 400 L. = l.KMAXl 
X = KMAX - L + 1 
X * DT / X 
DO 300 I = i.N 
DO 300 J = 1 » N 
PhH l.J» = 0.000 

c 

C PHI IS USED AS TEMPORARY STORAGE FOR THE CALCULATION OF GAMAT 

C 

DO POO K = 1 , N 
Y = F(I.K) * GAMAT ( K, J ) * X 
C 

C THE FOLLOWING FOUR STATEMENTS TEST FOR UNDERFLOW CONDITIONS. 

C 

IF (Y . fcQ. O.UDO) GO TO 2U0 
IF (ABS(Y) ,t,T» 1 .00—2 b ) GO TO 20 J 
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Y - G .01)0 
NN = NN + 1 

?J0 PHK l . J) » PHKI.JI + Y 

IF (I .£0. J) PHKI.J) = PHKI.J) + DT 
TOO CONTINUE 

on 40 J I = l.N 
on 4uo j = l.N 
40 C GAMATM.J) = PHKI.J) 

c 

L UPON CCMPLET ICN OF ThE L 00 LOOP, GAMAT CONTAINS THE DESIRED 

f. RESULT. THE REMAINDER OF THE SUBROUTINE COMPUTES PHI BY MULT I- 

C PLYING GAMAT bY F AND ADDING THE IDENTITY MATRIX. 

C 

nn 600 i * i.n 

DO 600 J = l.N 
PHKI.J) = O.UDO 
DO 600 K = l.N 

Y a F ( I . K ) * GAMATIK, J) 

IF (Y .ED. O.UDO) GO TO 6uj 
IF (AriSIY) .UT. 1.0U-2S) GO TO 600 

Y = 0.01)0 

• M N •= N N + 1 

YOU PHI (T ,J) = PHI ( l , J) + Y 

IF (I .FU. J ) PHKI.J) = Phi (l , J) + l.ODO 
bUC continuf 

IF I IOP .FO. 0) GO TO ttOU 

c 

r. PRIM output results 

c 

WRITE (6,700) NN 

700 FORMAT (IX / IX, 1 JHUNOE RFLU WS = . 1 5 // 19X . 3HPHI ) 

00 6GB J = 1,6 
6U8 FM(J) = VFMT(N.J) 

wRITF (o.FMT) ((PHKI.J), J = l.N), I = l.N) 
wRITF (6. 720) 

720 FORMAT (19X, 6HGAMAT ) 

DO 509 J = 1.6 
609 FMT(J) = VFMT(N.J) 

WRITE (b.FMT) ( ( GAMAT ( I , J) , J = l.N), I = l.N) 

BOO CONTINUE 
RETURN 
FND 
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itPFTC LUG LIST 

SUBRUUTJNF 01 AG {PHI, G AMI • uAMp, TO, PH 10 , VI, VP, N , C,N2, AAA, 

1 EIG, CPR , CPI. CR, LI, AH, LwVU, MWVU, RR, RI, PPHID, PHIP, TT, 

2 PXTlb, PERN, I PERN. IPER. OU , 0 , V, EXT2, SO, X, uS, CUU, PER2N, 

4 I PER P , IPEkPiM, PPS, PP, UP, USi, GF2, YD2 , IuP) 

C 

C 4444444444444444 4444444444444 44 4444 44 44 44 44 44 4444 44 4444444 444444444*44 

c 

C 1)1 Ati COMPUTES THE BLOCK U I AGONAL EURM (PHID) CF THE MATRIX PHI. 

C IT ALSO FINDS THE TRANSFORMED VERSIONS ( VL AND VP ) GF TWO 

L INPUT MATRICES ( GAM I AND GAMP RESPECTIVELY >. THE TRANSFORMATION 

f. MATRIX TO IS ALSO FOUND. 

C 

C 0 1 AG CALLS SUBROUTINES HSBG, EIGVEC , EIGUR, AND DMINV. 

C 

C 44 4 44 4 4444 444444 44 444 44 44444 4 44 4 44444 4 444 44 44 4 444 44444 4 44 4 44444444444 4 

C 

f DIMENSION OF INPUT VARIABLES LF LIST 

C PHi(N.N) 

i UAMi(N.C) 

C l, A M21 N.C ) 

C. O IMENS ION OF OUTPUT VARIABLES OF LIST 

C T C ( N , N ) 

C PE ID(N.N) 

C II N . c » 

C VP(N.C) 

C U I ME A S I UN OF INTERNAL VARIABLES 

C TT1 N , N ) 

f. RR 1 NJ 

r. run) 

C CR(N) 

C CUN) 

C AHIN.N) 

C LhVDlN) 

C MRVD(N) 

C PPhlD(N.N) 

C P H I P ( N . N ) 

C AAA(N.N) 

C CPRIN) 

r cp i ( n ) 

C F I G I N. N , ? > 

C FXT16lN.NI 

C OD(N.N) 

C 0 ( N , N ) 

C UUC(N.l) 

C OSIN.)) 

C SClN.l) 

C YIN, 1 ) 

C X ( N , 1 ) 

C PERNIN) 

c ipfknim 

C IPFR(N) 

C FXTPI2N.2N) 

C P2SI2N.2N) 

C P2I2N . PN) 

C DSl(PN.l) 

C 0P2I2N.I) 

C YD2I2N.1) 

r. opiPN.n 
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n r> n ^ r> r> r> n n 


C /2(2N,1) 

C PEK2NI 2N> 

C IPFR2NI2N) 

C (PFR2I2N) 

JNTEGFR C 

DOUBLE PRECISION PHI, GAM1, GAM2, TO, PHlD, VI, V2, PHIP, PPHID, 
1 IT. ObT 

0 IMFiviS ION PHIIN.N). GAMKN.C). GAM2IN.C) 

01 MENS ION TCIN.N). PHlD(N.N), V1IN.C). V2(N,C) 

0 1 WENS ION RR 1 N > » RI(N), GRIN), C1(M. AH(N.N), LWVD(N), MWVD(N), 

1 PPHiniN.NI, PnIPIn.N), AAA ( N ,N ) , CPRIN), CPIIN), EIGIN.N.2) , 

2 EXTJ6IN.N), UO(N.N), D(N.N> . OGQIN.I), OS(N.l), SOIN.l), YIN, I), 
■i X ( N » 1 ). PERNIN), I PERN ( N > , IPFR(N), EXT2IN2,N2), P2SIN2,N2), 

4 P2IN2.N2) . 0S1 I N2 » i ) » QP2IN2.1), YD2(N2,1)» U21N2.1) , 22IN2, 1) , 

5 PFR2NIN2) , IPER2NIN2). IPER2IN2), TTIN.N) 

COMMON / FORM / vFMTUO.o), FMT I 6 ) 


00 170 I ■= l.N 
CPR( I ) = 0.0 

cpim = o.o 
nr leg k = i,i. 
vni.Ki = o.ooo 

169 V2I I,K) = 0.000 
on 170 J = l.N 
PFIDI I .0) = 0.000 
AAAI I.J) = PHI 1 1 , J ) 

170 AH { I.J) = PHI! I.J) 

CALL HSHGIN, AH, N) 

OFTERMINE EIGENVALUES OF PHI MATRIX 

IF IIOP ,EQ. 0) GO TO 171 
WRITE (6.438) 

171 CALL EIGOk I AH, N, CR , Cl, lOP) 

00 290 II = l.N 
CPR I I I) = CR III ) 

CP HI I ) = CI< II ) 

290 CENT I NOE 

1 = 1 

5 PH 1 0.1 L » L 1 = CPR I L ) 

L = L + L 

IF (L .GT. N) GO TO 430 
IF ICPl(L-l) • EQ • 0.0) GO TO 5 
PHIOIL.L) = CPRlL-l) 

PHIOI L-l.L) = - CPIIL-l) 

PHlOIl .L-l) = CPI I L - 1 ) 

I = L + l 

IE (L .GT . ivl ) GO TO 9 30 
GO TO 5 

FORM EIGENVECTORS 

930 CALL EIGVEC (AAA, CPR, CPI, EIG. N, IOP, N2 , EXT16, PERN, IPEP-N, 

1 IPER, Dl». D , OS, 000, Y, SU. X, EXT2 , PER2N, I PE R2N , IPER2, P2S , 

2 P2 , CS1 , 02, UP2 , Y02, 12 ) 
no 935 I - l.N 

nr. 935 j = i,n 

TO (I.J) = ElGIl.J.l) + EIG { I.J, 2) 
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n r> r> 


435 


510 
44 5 


446 


4 50 

436 

437 
43 8 

439 

440 


AM I , J) = PH ID ( I « J) 

TTU.J) = TD( |,4) 

IF ( I0P ,EQ. 0) GO TO 4-»i 
-JRITF (6.436) 

on 5io j = 1,6 

FMKJ) * VFMT(N.J) 

»JkITF I6.FMTI ((PHIO(I.J), J = 1,N), I = 1,N> 

CAli HSriG (N, AH, N) 

DETERMINE FIGENVALUFS OF PHlO MATRIX 

I F ( IOP .EQ. 0) GO TO 446 
WRITF (6.439) 

0 ALL eiGOK (AH, N, HR • kl, IOP) 

CALL CM IN V (TT, N, OET , LWVU. MwVD) 

IF (OET .FO. O.OUO) WRITE (6,440) 

DO 450 K = l.C 
i)0 450 I = 1 *N 
on 450 J = 1 * N 

V1(I, K) - TT(I.J) * GAMKJ.K) + VUI.K) 
v/P ( I « K ) = TT(I.J) * GaM2(J.K) + V2( I *K) 

FORMAT ( 1 8 X , 4HPHID) 

FORMAT (IX. iP 10E 1 3 • 5 ) 

FDkMAT (IX. 23HE IGENVALUES UF PttI MATRIX) 

FORMAT (IX, 2oHE IGENVALUES OF PhID MATRIX) 

FORMAT (IX. 62FUET = 0.0 AND a LOOK DIAGONAL TR ANiFOKMAT I UN MATRIX 
1 IS SINGULAR ) 

RETURN 

FNU 


73 



r> n r> o n r> n n n r» n r> r> r> r> n n n 


SIBFTC EiGOR LIST 

SUHROUT INE EIGOR I a. M, ROOT k * ROOT 1 1 IPRNT) 

********************************************************************** 

FIGUR DETERMINES THE EIGENVALUES GF A REAL MATRIX <AJ USING 
FRANCIS OR ALGORITHM. MAXIMUM NUMBER OF ITERATIONS IN THE 
OR PROCESS IS SO. IF IPRNT IS NOT ZERO » THE EIGENVALUES 
WILL BE PRINTED OUT. 

********************************************************************** 

VARIABLES OF LIST 
AIM.M) INPUT MATRIX 

ROOTR I M« l > REAL PART OF EIGENVALUES 

ROOT l I M» I ) IMAGINARY PART OF EIGENVALUES 

DIMENSION AIM.M). ROOTR I M) » KOGTt(M), PS1I2), GRI3) 


N = M 

DO 2 I - l.N 
ROOT R I I ) = 0.0 
2 ROOT III) = 0.0 

IF < IPRNT) BO. 81 .BO 

80 WRITE 16.104) 

81 ZERO = 0.0 
OJ = I 

177 XNN = 0.0 
XN2 = 0.0 
A A = 0.0 
B = 0.0 
C = 0.0 
00 = 0.0 
R = 0.0 
SIC- = 0.0 
ITER = 0 

17 IF IN- 2) 13,14.12 

13 IF I IPRNT) 82,83,82 

82 WRITE (6,108) A(l.l) 

8 3 ROOT R I 1 ) = A (1.1) 

ROOT I ( 1) = 0.0 
1 RETURN 

14 JJ = - l 

12 X = (A(N-l.N-l) - A IN ,N) ) ** 2 
S = 4.0 * AIN.N-1) * AIN-l.N) 

ITER = ITER + l 

IF (X .FO. 0.0 .UK. ABSIS / X) .GT. l.OE-8) GO TO 15 

16 It- I ABSIAI N— l.N— 1) ) - ABSIA1N.N) ) ) 32,32.31 

31 F = AIN-l.N-l) 

G = A(N.N) 

GO TO 33 

32 G = A 1 N— 1 , N— 1 ) 

F = AIN.N) 

33 F - 0.0 
H = 0.0 
GO TO 24 

15 S = X + S 
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X = AIN— i»N— 1) + A(N»N) 

I F IS) 18,19,19 

19 SO = SORT(S) 

F = 0.0 

H = 0.0 

IF (X) 21,2 1,22 

21 F = (X - SO) / 2.0 

r, = (X ♦ SO) / 2.0 

Gfl TO 24 

22 G = (X - SO) / 2.0 

F = (X + SO) / 2.0 

Gfl TO 24 

18 F » SORT ( - S) / 2.0 
F = X / 2.0 
li = E 
H = - F 

24 IF CJJ) 2 8,70,70 

70 i) = 1.0E-10 * (ABS(G) + F) 

IF (ABSIAI N-l.N-2)) .GT. 0) GO TO 26 

28 IF (IPRNT) 84,85.84 

84 WRITE (6,105) E, F, ITER 
WRITE (6,105) G, H 

85 ROOTRIN) = E 
kOOTUN) = F 
ROOTR(N-i) - G 
ROOT I ( N- 1 > * H 
N = N — 2 

IF (JJ) 1,177,177 

26 IF < ABSI AIN, N— 1 ) ) .GT. 1.0E-10 * ABSI A( N, N) ) ) GO TO 50 

29 IF ( IPRNT) 86,87,86 

86 WRITE (6,105) A ( N , N ) , ZERO, ITER 

67 ROOTRIN) = AIN.N) 

ROOT I (N) =0.0 
N = N - 1 
Gfl Tfl 177 

50 IF ( ABS( AB S( XNN / A(NtN-l))- 1.0) - i.GE-6) 63.68,62 

62 IF ( ABSI AtiS( XN2 / A(N-i,N-2)) - l.U) - 1.0E-6) 63,63,700 

63 VO = ABS(A(N, N— 1 ) ) - ABS ( A< N- 1 , N-2 ) ) 

IF ( ITFR - 15) 53,164.04 

1 64 IF (VO) 165,165,166 

165 R = AIN-1.IM-2) 2 

SIG - 2.0 ♦ A( N-l.N-2) 

GO Tfl 60 

166 R = AIN, N— 1 ) *4 2 
SIG = 2.0 4 AIN.N-l) 

GO Tfl 60 

64 IF (VO) 67,67.66 

66 IF I IPRNT) 88, 85,88 

88 WRITE 16,107) AIN-l.N-2) 

GO TO 84 

67 IF I IPRNT) 89,87.89 

89 WRITE (6,107) AIN.N— 1) 

GO TO 86 

700 IF (ITER .GT. 50) GO TO 63 
IF I ITER .GT. 5) GO TO 53 

701 /. 1 = (It - AA) 4 * 2 + (F - B) 44 2) / IE 4 £ + F 4 F) 

72 = (10 - C) 4* 2 + (H - DO) 4* 2) / (G*G + H4H) 

IF ( Z 1 - .25) 51.51.52 

51 IF (72 - .25) 53.53,54 

53 R = E*G— F*H 
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SIG = E + G 
GO TO 60 

54 R = E * E 
SIG = E + E 
GO TO 60 

3^ IF 1/2 - .25) 56,55,601 

55 R = G * G 
SIG = G + G 
GO TO 60 

001 R a 0.0 
SIG = 0.0 

6C XfcN a AIN.N-1) 

XiM2 = A ( N— 1 » N— 2 ) 

N1 * N — L 
IA * N - 2 
IP * IA 

IF IN - 3) 201 .2 10 » 260 

260 00 212 J = 3*N 1 
,11 - N — J 

IF (AfiSlAI Ol + i, J1)J - 0) 210.210,211 

211 OEN = AlOk+l.Ol+l) * (A<01+1,01+1> - SIG) + AIOl+1,01+2) * 
l AI01+2.01+1) + R 
IF I OFN ) 261,212.261 

261 IF ( AHStAI Jl+1, 01) * AI01 +2,01+1) * ( AB S( A( 01 +1 , J 1+1 ) + 

l A(0 1+2,01+2) - SIG) + ABSIA101+3.01+2) >> / OEN) - 0) 210,210,212 

2 12 IP * 01 

210 on 214 0 = l.IP 

01 = IP - J + 1 

IF 1 ABSIAI 01+1,01 )) - 0) 213,213.214 

214 10 a J1 

213 00 20C I = IP.Nl 

IF I I - IP) 216.215,216 

215 GR(l) - A I IP, IP) * I A I IP, IP ) - SIG) + A(IP,IP + i) * At IP+1, IP) + R 
GR I 2 ) = A ( IP + 1, IP ) * <A< IP, IP) + A( IP + 1, IP+1) - SIG) 

.;R(3> = A( IP + 1, IP) * At IP +2, IP+1) 

A ( IP +2, IP) = 0.0 
GG TO 215 

2 16 GRI J ) - A( I, 1-1) 

GR(2 ) - A ( 1+ 1, 1-1 ) 

IF I l - IA) 217,217,216 
217 GR ( 3 ) = A ( 1+2. l-l ) 

GO TO 215 
21b GRI 3) = 0.0 

215 XK = SIGN! SURTlGR(l) ** 2 + GRI2) ** 2 + GR ( 3 ) ** 2), GR(ll) 

222 IF IXK) 2/3,224.223 

223 At = GRI 1 ) / XK + 1.0 
PS1U) = GRI 2) / (GR ( 1 ) + XK) 

P S 1 1 2 ) = GR I 3 ) / (GR(l) + XK) 

GO TO 225 

224 At = 2.0 
PSIM ) * 0.0 
PS 1 1 2) = 0.0 

225 IF I t - 10) 226,227.226 

226 IF I I - IP) 225,228,22V 
228 A I I. l-l) = - A( 1,1-1) 

GO TO 227 

225 A I 1,1—1) = - XK 

22/ 00 230 0 - I , N 

IF I I - I A ) 231,231, 232 

231 CR = P S I I 2 ) * A (I +2,0) 
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23 2 
2 33 


234 

230 

235 

2 36 

231 

238 

239 
241 


242 

240 

24 3 


20 G 
2 01 


104 

105 
10 7 


GO TO 23 3 
OK = 0.0 

FR = AL * (A( 1* J) + PSi(i) * AU + i.J) + CR) 

AUtJ) = AUtJ) ~ ER 

AU+UJ) A(lHtJ) - PSI(i) * ER 

IF { I - I A > 234 . 234 « 230 

A(I+;.J) = AU+2,J) - PSI (2) * ER 

CONTINUE 

IE i l - I A J 235*235,236 
L = I + 2 
GO TO 231 


L - N 

on 240 J - IQ.L 

IE ( i - I A ) 236 , 236*239 

CR - PSi(2i * AIJ.I+2) 

GO TO 241 
CR - U.U 

ER = Al * (A(J*I) + PS l ( 1 ) * A ( J • I + 1 ) + CR) 

A { J * I ) ~ A { J « l ) - ER 

AlJ.I+1) A(J.l+l> - PSI M> * ER 

IF <1 - I A) 242 • 242 * 240 

A ( j « I +2 ) = A (J* 1 + 2) - PSI ( 2 ) * ER 

CONTINUE 

IF (I - N +3) 243.243.200 

ER * AL * PS i ( 2 ) * A { 1+3 * 1+2 ) 

A I 1+3. I ) = - ER 

a 1 1 + 3. i + i> ~ - psu i> * er 

A i 1+3 * l + 2 ) = A( 1 + 3. 1+2) - PSI ( 2 ) * ER 
CENT INUE 
A A = E 
R = F 
C = G 
00 = H 


GO TO 12 

FORMAT {////IX. 9HREAL PART. 6X. I4H I MAG INARY PART. 
1 UhTAkEN AS ZERO. 6X. 4HITER //) 

FORMAT (IX. El 5. 8* 3X. E13.3. 42X. 13) 


FORMAT { 56X . E13.b) 
END 


26X, 
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SI8FTC F IGV EC LIST 

SUBROUTINE EIGVEC (AAA, CPR . CPI. El G. N, IGP2. N2. EXT 16* PERN, 

1 I PERN * I PER * DU. U* GS. GQG • Y, SO* X* EXT2 * PER2N* i PER2N* 

2 I PER 2 * P2S * P2 , GS1 * G2 * 0P2 * YD2* Z2I 
C 

C 

C F I GV EC DETERMINES TEE EIGENVECTORS OF A REAL MATRIX (AAA), 

C GIVEN THE REAL AND IMAGINARY PARTS OF THE EIGENVALUES OF (AAA), 

C LPH AND CPI. RESPECTIVELY. THE TECHNIUUE USED IS THE VAN NESS 

C INVERSE ITERATION METHOD. EIGENVECTORS ARE STORED IN A TRIPLE 

C SUBSCRIPTED ARRAY { £ I G) . THE REAL PARTS ARE STORED COLUMNWISE 

C IN EIGU.J.l). AND THE IMAGINARY PARTS IN EIGII.J.2). IF I0P2 

C IS NOT ZERO. THE EIGENVECTORS ARE PRINTED OUT. 

C 

C EIGVFC CALLS SUBROUTINES FACTR. PERM. AND GMPRD. 

C 

£ *#$$#&###&&#******#************#***#*#*#*#**#*#**##*#*#**#************ 


C 

c 


c 

C 

c 

c 

c 

c 

C 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


DIMENSION OF INPUT VARIABLES CF LIST 
AAA(N.N) 

CPR ( N ) 

CP I(N > 

0 1 MEN S ION OF OUTPUT VARIABLE CF LIST 
FIG (N.N.2) 

DIMENSION OF INTERNAL VARIABLES 
EXT 16 I N. N) 

PERN ( N ) 

IPERNI M 
IPER(N) 

DGIN.NI 
OIN.N I 
OSIN* 1 ) 

UQUIN * 1 ) 

Y IN* I I 
SO IN* I ) 

X IN* i I 
EXT2I 2N.2N) 

P2SI 2N.2N) 

P2I2N* 2N) 

QSlIPN.l) 

UP 2 1 2N * 1 ) 

Y02I2N. 1) 

Q2I2N* 1) 

/ 2 I 2N . I ) 

PFR2NI2N) 

IPER2M2N) 

IPFK2I 2N) 


DIMENSION AAAIN.N) . CPRIN). CPHN) 

D [MENS ION E IGIN.N.2I 

DIMENSION EXT 16 1 N. N) . UOIN.NJ. DIN, Nit YIN.ll, SDIN.i), 

L X ( N * 1 ) * YD2 ( N2 * 1 > » Z2(N2.l>. P2(N2,N2>, 0Sl(N2,i>, 0P2(N2»1> 

2 PERN(N>, I PERN (N ) • IPERIN). OS(N.l>, QUHN.IJ » Q2(N2,i). 

3 EXT2(N2«N?) * PER2NCN2), IPER2N(N2), IPER2(N2). P2SCN2.N2) 


? 
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L = 1 

310 IF CL .GT. N > GO TO 430 

IF (CPICL) .NE. 0.0) GO TO 360 
C CALCULATE REAL EIGENVECTOR 

DO 313 I = L.N 
DO 315 J = l.N 
FXT16C I, J) = AAA ( I. J) 

IF Cl .El). J) FXT16CI.J) = EXT Lb C 1 » J ) - CPRC L) 

315 CONTINUE 

316 CALL FAC TR ( EXT16 * PERN, N , N , IER) 

DU 313 I = l.N 

313 l PERN C I ) = PERN ( I ) 

CALL PERM CIPERN, IPEK, hi , IERP) 

IF ( IERP .EO. 0) GO TO 312 
*RITE (6,311 ) 

311 FORMAT ( 1 HO , 19HPFRN CANNOT BE DONE) 

312 CCNTINUF 

IF (IER .NF. 3) GO TO 315 
ORITE (6.317) 

317 FORMAT C1H0, 14HFACTK IS WRUNG) 

DO 318 l = L.N 

DO 31b J = l.N 

IF (EXT16CI.J) .Eli. 0.0) EXT16CI.J) = .5 / C2.0 ** 33) 

318 CONTINUE 
GO TO 316 

319 DO 320 I = l.N 
on 320 J = l.N 
DDC I , J ) = 0.0 
DC I, J) = 0.0 

IF C I .EO. J) DDC I , J ) = 1.0 

320 CONTINUE 

DO 323 I = l.N 
I I PER = I PERI I) 

DO 325 J * l.N 
DC I.J) a ODC IIPFR.J) 

325 CONTINUE 

DO 340 I = l.N 

IE (EXT 161 I « I ) .NE. 0.0) GO TO 340 
AMAX = Art SC E XT lb C I, 1) ) 

DO 335 J = l.N 

IE CA8SC EXTlbC l, J ) ) .GT. AMAX) AMAX = ABSCEXT16C I.J) ) 
335 CONTINUE 

IE (AMAX .EO. 0.0) AMAX = 1,C 
FXT16CI.C) = .5 * AMAX / 12.0 ** 35) 

340 CONTINUE 

DO 345 I = l.N 
OSC I.J) a l.U 

345 OOUI 1,1) = OSC I .1 ) 

.JO 1 END = 0 

346 CALL GMPRU (0, 000, Y, N , N , 1) 

SDCl.l) = Y(l.l) 

DO 34 7 I = 2.N 
SDCl.l) a Y( 1.1 > 

M = I - 1 
DO 347 J = l.M 

34 7 SDCl.l) a so l 1,1) - EXTlbC I.J) * SOCJ.l) 

X ( N .1) = SDCN .1) / EXTlbC N ,N ) 

00 349 I = 2 » N 
J = N - I + 1 
XCJ.l) a SDCJ.l) 
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M = J + 1 
on 348 K = «,N 

a48 X(J,1) = X(J,1) - X ( K , 1 ) * EXT16(J,K> 

349 X ( J . 1 ) - XtJ.l) / EXT16IJ.J) 

1/ - 1 

/MAX = A8S(X(1,1)) 

00 350 I = 2.N 
7 INT = AbS(X(I, 1) > 

IF ( Z I NT .IE. /MAX) GO TO 350 
/ MAX = Z INT 
1/ = I 

350 CONTINUE 

/MAX - l .0 / X( IZ.l) 

00 351 l = 1,N 

351 000(1,1) = X(I.l) * ZMAX 
10 S o 

DO 35 5 l = UN 

IE 1 ArtS (000( 1*1)) .GT. UOE-ld) GU TO 352 

0001 I. 1) = 0.0 
GO TO 353 

352 0 I E = AdS( uS ( l , l ) - 000(1 ,1)) 

IF (DIF .GT. .001 * AbS ( 000 (1,1)) ) GG TO 355 

353 10 = 10 + i 

355 CONTINUE 

IF (10 .EO. N ) GO TO 420 
J01END = JO 1 END + 1 
IE ( JU1END .EO. 20 GO TO 358 
DO 356 I = l.N 

356 US(I.l) * 0001 1.1) 

IF (J01FND .LT. 30) GO TO :>46 
wRITE(e.:>57) CPR(L) 

357 FORMAT (IX, 29HSTUCK IN LOOP WHERE CPR(L) = , 

12oEAnD ANSWERS MAT NUT BE RIGHT) 

GO TO 420 

358 00 359 I * l.N 

uoo(l.i) = (000(1.1) + OS ( I , 1 ) ) / 2.0 

359 O S ( I , 1 ) = COO ( I . 1 ) 

GO TO 346 

C CALCULATE COMPLEX EIGENVECTOR 

360 DO 354 I = l.N 

DO 354 J = l.N 

354 FXT 1 6 ( I.J) = AAAI l.J) 

DO 364 I = l,N2 

DO 364 J = i,N2 

564 EXT2 ( I.J) = U.O 
DO 36 5 I = UN 

DO 365 J = l.N 

K = I + N 
J J •= J + N 

FaT2 ( I.J) = EXT 16(1 . J) 


FXT 

,JJ ) 

■= 

EXT 16 ( I , J) 



IF 

i 1 

• FU • 

J > 

EXT2 ( K , J ) 

= - CP I ( L ) 


IF 

( i 

• E Q * 

J i 

EX T2 ( I .JJ) 

= CPI ( L ) 


l F 

C l 

*E u) * 

J ) 

EXT2 ( I.J) 

= EXT16( I.J) 

- CPR(L) 

If 

< I 

• FiJ • 

J > 

EX T2 ( K , J J ) 

= E XT 16 ( I ,J) 

- CPR( L) 


365 CONTINUE 

366 CALL EACTR ( EXT2 , PEK2N, N2 . N2 . IER) 
DO 363 l - UN 2 

363 I PER 2N ( i ) = PER2N ( I ) 

CALL PERM (IPER2N, IPERz. N2. IERP) 


1PE14.6, 
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IF < l ERP .EQ. 0) GO TO 362 
WRITE (6.361) 

361 FORMAT (1H0, 20HPER2N CANNOT BE DONE) 

362 CONTINUE 

IF ( IER .NE. 3) GO TO 369 
WHITE (6.367) 

367 FORMAT (lHO, 16HEACTR 2 IS WRONG) 

00 368 I = l,N2 

00 368 J ■= 1.N2 

IF (FXT2( t.J) .EU. 0.0) EXT2U.J) = .5 / (2.0 ** 35) 

368 CONTINUE 

GO TO 366 

365 00 .370 I = 1.N2 

on 370 J = 1.N2 

P2S( I . J) =0.0 

P 2 ( I » 0 ) = 0.0 

IF (I .EU. J) P2SII.J) = 1.0 
370 CONTINUE 

on 375 I = 1 . N2 

1 I PER = IPER2( I ) 

Oil 37 5 J = 1.N2 
P21I.J) - P2 S( I IP ER . J ) 

375 CONTINUE 

no 390 I = 1.N2 

IF ( EXT2 ( I.I ) .NE. 0.0) GO TO 390 
AMAX = Art S ( E XT2 ( l .1) ) 

1)0 38 5 J = I.N2 

IE IABSI£XT2( I. J) ) .GT. AMAX) AMAX = A8S ( EXT2 ( I . J ) ) 
385 CONTINUE 

IF (AMAX .FQ. 0.0) AMAX - l.C 
EXT 2. (I.I) = .5 * AMAX / (2.0 ** 35) 

390 CONTINUE 

00 395 I = 1.N2 
US 1(1.1) = 1.0 

395 021 1,1) = CS 1(1,1) 

J02EN0 - 0 

396 CALL GMPRO (P2, 02, 0P2. N2 , N2, 1) 

Y02 ( 1 » 1 ) = 0P2U.1) 

00 397 I = 2.N2 
YQ2 ( 1.1) = UP2i I , i) 

M = I - 1 
On 397 J - l.M 

397 Y 02 ( 1.1) = YD2M.1) - EXT2U.J) * YD2(J,1) 

77(N2.1) = Y02(N2,1) 7 EXT2(N2.N2) 

00 399 I = 2.N2 
J = N2 - I + 1 
/ 2 ( J , 1 ) = Y I) 2 ( J , 1) 

M = J + 1 
on 398 K = M.N2 

398 7 2 ( J . 1 ) = Z 2 ( J , 1 ) - 72IK.1) * ExT2(J,K) 

395 7. 2 ( J » 1 ) * 72 ( J « 1 ) / EXT2(J,J) 

IMAX = 1 
J = N + 1 

/MAX = 72(1,1) * ip (l.l) + Z 2 ( j , 1 ) v 72 ( J ,1 ) 
on 400 I = 2.N 
J = I + N 

7 I NT = 72( 1. 1) * 72 ( 1 .1) + 72(J,1) * Z2(J,1) 

IF (ZINT .LE. ZMAX) GO TO 900 
/MAX = 7 I NT 
IMAX = 1 



400 CC.NTlN.ue 

UMAX = I MAX + N 
/MAX = SORT! ZMAX) 

00 401 I ■= l.N 

.) - I + N 

02(1,1) = (Z2U.1) / /MAX) * ( Z2 ( I MAX » 1 ) / ZMAX) + 1Z21J.1) / 

1 /MAX) * (Z21JMAX.1) / /MAX) 

02 ( J » 1 ) = - IZUI.il / ZMAX) * ( Z2 ( JMAX » 1 ) / ZMAX) + (Z2(J,l) / 
1 /MAX) * ( Z2 ( IMAX .1 ) / ZMAX) 

401 CL NT I NUE 
10 = 0 

OG 405 I = 1.N2 

IF ( A8S1 021 1*1) ) .GT. l.OE-lO) GU TU 402 
02(1.1) = 0.0 
GO TU 40.3 

402 OIF = ABSl USi( 1.1) - 02(1.1)) 

IF (OIF .GT. .001 * AbS (02(1.1))) GO TC 405 

403 10 = IQ + 1 

405 CCNTINUF 

IF ( 10 .EU. N2 ) GU TO 410 

.J02EN0 - J02EN0 + 1 

IF ( J02EN0 .Eu. 20) GO TU 40b 

00 4 C o I = 1.N2 

406 OSK I , 1 ) = 02(1.1) 

IF (.J02END .LT. 30) GO TU 396 
wRITE(6,40/) CPK ( U ) » CPI(L) 

407 FORMAT (IX. 29HSTUCK IN LUUF WHERE CPR(L) = , 1PE14.6, 

1 13HAMC CP I ( L ) = , iPE14. 6, 

228FANC ANSWERS MAY NOT ttfc RIGHT) 

GO TO 410 

408 00 40 9 I = 1.N2 

02(1,1) = (02(1.1) + CSl(t.l)) / 2.0 

409 OSK I . 1) ■= 02(1.1) 

GU TO 396 

410 CONTINUE 

C LOAD EIGENVECTOR MATRIX WITH COMPLEX EIGENVECTOR PAIR 

on 413 I = l.N 
.1 = l + N 

E I G ( I , L . 1 ) = 02(1.1) 

FIG( ( , L * 1 * 1 ) = 02(1,1) 

F IG( I , L . 2 ) = 02 ( J . 1 ) 

413 F IG( I , L+ 1 » 2 ) = - 02 ( J , 1 ) 

L = L + 2 
fin TO 310 

C END OF COMPLEX EIGENVECTOR CALCULATION 

C LOAD EIGENVECTOR MATRIX WITH REAL EIGENVECTOR 

420 00 425 I = l.N 

FIG(I.L.l) = OuOI i . I) 

E l G ( I ,L . 2) = 0.0 
425 CCNTINUF 

1 = L + 1 
fill TO 310 

C PRINT OUT EIGENVECTOR MATRIX 

430 LL = l 
LLL = N 

IF (N .GT. 10) LLL - 10 
IF ( I0P2 .FU. 0) GO TO 467 
WRITE (6.440) 

440 FORMAT (1H1 / 20X, 12hE l GENVEC TORS //) 

435 00 455 I -= l.N 
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WRITE (6,445) { t IG { l ,L , 1 ) » L ••» LL,LLL) 

445 FORMAT (JLX, IP10E12.4) 

WRITE (6,450) (E IGH ,L *^) , L * LL.tLL) 
450 FORMAT ( 5X * lP10Ei*>.4> 

455 CONTINUE 

WRITE (6,465) 

465 FERMAT ( LH1 ) 

467 TF (ILL .EC, N ) GO TO 460 
LL « II 
LlL = N 
GO TO 435 
460 CONTINUE 
RETURN 
FNO 


if BFTC PERM l 1ST 

SUBROUTINE PER M IP I , IP2, N, IERI 

****** ** 4c * * *** * ** * * *** *** * 4: 

COMPUTE PERMUTATION VECTOR IP*' FOR TRAN SP C SI T I CN VECTOR IPi 
******************************************* *************************** 
0 I MENS ION IPJLID* IP^(i) 


00 2 I s 1,N 
/ l P2 ( l ) - I 

on 6 I = i.n 
k * ipim 
IF ( K - I > 3.6,4 

3 IF (K) 7,7.3 

4 IF IN - K) 7,3,6 

5 .1 •= IP2I I ) 

IP2< I > = IP2IK) 

IP2IK) = J 

6 CONTINUE 
IFR « 0 
RETURN 

ERROR RETURN - IPi IS NUT A TRANSPOSITION VECTOR 

7 IFR * l 
RETURN 
END 
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